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This  research  focuses  on  the  control  of  a  nonlinear  system  whose  output  sub- 
ject to  an  additive  disturbance.  The  main  interest  is  in  the  investigation  of  controllers 
that  reduce  the  effect  of  the  disturbance  on  the  system  output.  Usually,  it  is  not  possi- 
ble to  construct  a  controller  that  completely  eliminates  the  effects  of  the  disturbance. 
It  is  then  of  interest  to  find  how  well  the  "best"  controller  can  attenuate  the  effect  of 
the  disturbance.  The  main  result  of  this  dissertation  is  a  performance  bound,  that 
provides  an  estimate  of  the  best  disturbance  attenuation  that  can  be  achieved  for  a 
given  system,  using  a  causal  controller  that  renders  the  system  internally  stable. 

An  approximate  right  inverse  of  a  nonlinear  system  E  is  introduced  to  fa- 
cilitate the  derivation  of  the  performance  bound  and  the  development  of  nonlinear 
controllers.  An  approximate  right  inverse  is  constructed  to  be  stable  and  causal  for 
implementation  purposes.  The  difference  between  an  approximate  right  inverse  and 
a  right  inverse  is  that  the  right  inverse  may  not  be  both  stable  and  causal.  The  role 
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of  an  approximate  right  inverse  is  to  approximate  a  disturbed  signal  with  a  signal  in 
the  image  of  the  system  E.  An  approximate  right  inverse  can  be  constructed  for  any 
system  E. 

The  calculation  of  the  performance  bound  involves  an  optimization  process  of 
finding  a  global  maximum  of  a  non-convex  function.  For  several  cases,  a  nonlinear 
programming  algorithm  is  developed  to  handle  the  optimization.  To  demonstrate  the 
application  of  this  performance  bound,  it  is  calculated  for  three  practical  systems: 

1.  the  voltage  control  of  a  permanent  magnet  stepper  motor; 

2.  the  longitudinal  control  of  an  aircraft;  and 

3.  the  multivariable  process  control  of  a  regulator  that  regulates  the  liquid  level 
in  a  pressurized  tank. 
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CHAPTER  1 
INTRODUCTION 


Over  the  last  2  decades,  there  has  been  considerable  interests  in  the  literature 
in  the  derivation  of  optimal  controllers  that  minimizes  the  effect  of  the  disturbance  on 
the  output  of  a  control  system.  This  dissertation  addresses  this  question  for  the  case 
of  nonlinear  systems.  The  main  objective  is  to  derive  an  estimate  of  the  performance 
of  an  optimal  controller,  by  deriving  a  bound  on  the  effect  of  the  disturbance  has  on 
the  system  output  when  the  optimal  controller  is  used.  Using  this  bound,  we  can 
then  gauge  the  performance  of  suboptimal  controller,  to  see  how  well  they  compare 
to  the  optimal  ones.  Suboptimal  controllers  may  be  much  easier  to  implement  than 
their  optimal  counterparts.  Thus,  our  performance  bound  can  be  use  to  find  simple 
controllers  whose  disturbance  attenuation  properties  are  close  to  those  of  optimal 
disturbance  attenuating  controllers. 

The  performance  bound  derived  in  this  dissertation  came  from  the  require- 
ment to  characterize  the  performance  of  a  controller  that  reduces  the  effect  of  the 
disturbance  on  the  output  to  a  minimum.  The  original  work  provided  is  the  derivation 
of  the  performance  bound  and  nonlinear  programming  algorithm  in  how  to  calculate 
the  performance  bound.  The  application  of  the  performance  bound  for  disturbance 
attenuation  is  proposed  as  future  directions. 

The  basic  design  problem  to  which  this  performance  bound  relates  is  the 
problem  of  disturbance  attenuation  for  nonlinear  control  systems.  Specifically,  this 
dissertation  discusses  the  following  configuration.  In  the  configuration  of  Figure  1.1, 
S  is  the  nonlinear  system  to  be  controlled.  C  represents  an  equivalent  controller  that 
incorporates  all  the  control  elements  of  the  loop.   The  external  (or  reference)  signal 


is  denoted  by  v;  the  disturbance  signal  is  denoted  by  d\  and  the  output  signal  is 
denoted  by  z.  The  closed  loop  system  is  required  to  be  internally  stable.  Internal 
stability  signifies  that  a  configuration  can  tolerate  small  disturbances  on  its  external 
and  internal  ports  (including  ports  within  the  equivalent  controller  C)  without  losing 
stability.  The  equations  that  describe  Figure  1.1  are 


z  =  d  -f  y,      y  =  Eu,      u  =  C(v,z). 


(1.1) 


In  Figure  1.1,  He  represents  the  appropriate  equivalent  system.  This  can  be 
expressed  in  notation 

z  =  Y,c{v,d)  (1.2) 

where  the  output  signal  z  is  determined  by  the  signals  v  and  d  and  depends  on  the 
system  E  as  well  as  the  equivalent  controller  C .  We  would  like  to  reduce  as  much 
as  possible  the  effects  of  d  on  z.  Our  bound,  which  is  the  attainable  performance 
for  control  of  a  nonlinear  system  whose  output  is  subject  to  an  additive  disturbance, 
provides  an  estimate  of  the  minimal  effect  of  d  on  z.  Using  the  estimates  of  the  min- 
imal effect,  we  can  evaluate  controllers.  The  analysis  of  the  dissertation  is  restricted 
to  the  case  of  discrete-time  systems. 

The  desired  response,  with  or  without  the  disturbance  signal  for  the  configu- 
ration in  Figure  1.1,  is  z  =  Eu.  To  null  out  the  disturbance  signal,  the  design  of  the 
equivalent  controller  C  would  be  such  that  y  =  Eu  =  Eu  —  d.   Substituting  for  y  in 
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Figure  1.1:  The  block  diagram  of  Ec. 


the  first  equation  of  (1.1),  this  yields  z  —  Ei>  which  is  our  stated  goal.  For  the  sake  of 
argument,  assume  that  the  system  E  has  a  stable  and  causal  inverse  system  E_1.  The 
input  signal  u  to  the  system  E  would  then  be  required  to  be  u  —  E_1(Eu  —  d).  The 
inverse  system  E_1  in  some  cases  might  not  be  implementable  because  it  is  not  sta- 
ble or  it  is  not  causal.  For  those  cases,  approximations  of  the  inverse  systems  would 
then  be  used  for  practical  implementation.  The  approximations  will  be  selected  in  a 
manner  suited  to  approximate  the  signal  T,v  —  d  with  a  signal  in  the  image  of  E. 

The  composition  of  the  system  and  its  approximate  inverse  system  forms  a 
nearly  identity  system.  For  points  inside  the  image  of  the  system  E,  the  composition 
will  appear  as  the  identity.  For  external  points,  which  are  points  outside  the  image  of 
the  system  E,  the  composition  will  produce  a  point  in  the  image  of  E  that  is  closest 
to  the  external  point.  Our  performance  bound  is  a  norm  that  gauges  the  closeness  of 
the  identity  system  and  a  nearly  identity  system  formed  from  the  system  S  and  its 
approximate  inverse. 

The  performance  where  E<?  can  attain  by  using  our  bound  is  determined  by  the 
inherent  properties  of  the  system  E.  The  bound  is  directly  related  to  an  approximate 
inverse.  The  role  of  an  approximate  inverse  is  to  approximate  the  signal  u  with  a 
signal  v  in  the  image  of  E.  Therefore,  the  approximate  invertibility  of  the  system 
E  provides  a  measure  of  the  ability  to  match  a  set  S  by  a  subset  of  the  image  of 
the  system  S.  A  measure  of  singularity  is  introduced  as  an  indicator  of  approximate 
invertibility  which  is  then  used  to  derive  the  performance  bound.  To  demonstrate  the 
application  of  this  performance  bound,  it  is  calculated  for  three  practical  systems: 

1.  the  voltage  control  of  a  permanent  magnet  stepper  motor; 

2.  the  longitudinal  control  of  an  aircraft;  and 

3.  the  multi variable  process  control  of  a  regulator  that  regulates  the  liquid  level 
in  a  pressurized  tank. 


The  dissertation  is  organized  as  follows.  Chapter  2  contains  the  basic  notions 
of  nonlinear  systems  including  the  theory  of  fraction  representation.  Chapter  3  dis- 
cusses an  approximate  right  inverse  and  the  performance  bound.  Chapter  4  discusses 
the  main  results  of  the  dissertation,  including  calculations  of  the  performance  bound 
for  applications.  Chapter  5  contains  a  summary  and  future  directions. 

1.1     Background 

This  section  contains  a  qualitative  survey  of  some  topics  in  nonlinear  control 
that  are  important  to  the  dissertation.  A  more  technical  discussion  of  these  topics  is 
provided  in  Chapter  2. 

The  description  of  £c  in  Figure  1.1  has  its  origins  in  the  theory  of  fraction  rep- 
resentation of  nonlinear  systems.  It  can  be  stated  that  a  right  fraction  representation 
of  a  nonlinear  system  £  is  a  factorization  of  £  into  a  composition  form  E  =  PQ~l , 
where  P  and  Q  are  stable  systems  with  Q  being  invertible.  In  Hammer  [16,  18],  tools 
for  a  compact  and  perceptive  statement  of  results  were  developed  for  the  theory  of 
fraction  representation  of  nonlinear  systems.  The  fraction  representation  £  =  PQ~l 
is  said  to  be  coprime  when  the  systems  P  and  Q  are  right  coprime.  An  attribute  of 
right  coprime  fraction  representation  is  that  every  instability  of  the  inverse  system 
Q~l  is  also  an  instability  of  the  system  E;  i.e.,  there  is  no  cancellation  of  instabilities 
within  the  composition  PQ~l . 

A  causal  (respectively,  strictly  causal)  system  £  is  one  where  the  values  of  the 
output  sequence  Eu  up  to  and  including  index  i  (respectively,  i  -f  1)  depend  only  on 
the  values  of  the  input  sequence  u  up  to  index  i.  A  system  is  bicausal  if  it  is  causal 
and  if  it  possesses  a  causal  inverse. 

In  this  dissertation,  it  is  assumed  that  the  system  E  being  controlled  can  be 
stabilized,  that  it  is  strictly  causal,  and  that  it  possesses  a  right  coprime  fraction 
representation  of  the  form  E  =  PQ~l,  with  Q  being  bicausal.     The  stabilization 


assumption  on  the  system  £  is  necessary  because  the  closed  loop  system  is  required  to 
be  internally  stable.  Strict  causality  is  placed  on  the  system  E  as  a  handy  assumption 
to  certify  that  the  closed  loop  system  is  well  posed.  Strict  causality  is  not  an  essential 
condition  and  it  can  be  replaced  with  plain  causality  combined  with  a  well-posedeness 
requirement.  This  dissertation  relies  on  stabilization  theory  which  applies  only  to 
systems  possessing  right  coprime  fraction  representations. 

A  result  from  [22]  gives  a  simple  parameterization  of  the  set  of  all  system 
responses  that  can  be  obtained  through  internally  stable  control  of  a  given  system. 
The  control  scheme  to  control  E  is  shown  Figure  1.1.  The  parameterization  provides 
a  clear  indication  of  the  effects  of  the  disturbance  on  the  response  of  the  stabilized 
closed  loop  system.  This  result  can  be  summarized  as  follows.  Let  E  =  PQ~l  be 
a  right  coprime  fraction  representation  (with  a  bicausal  "denominator"  Q)  of  the 
system  being  controlled.  Then, 

1.  For  every  causal  equivalent  controller  C  for  which  the  closed  loop  system  of 
Figure  1.1  is  internally  stable,  there  exists  a  stable  and  causal  system  <f)(v,d) 
such  that 

Y,c(v,d)  =  d  +  P(f>(v,d)  =  [I  +  P<j>(v,-)]d  (1.3) 

where  P  is  the  "numerator"  of  the  right  coprime  fraction  representation  of  £ 
and  /  denotes  the  identity  system. 

2.  Conversely,  for  every  stable  and  causal  system  (f>{v,  d),  there  is  an  internally  sta- 
ble control  configuration  around  the  system  £  for  which  the  equivalent  system 
£c(u,  d)  satisfies  £c(u,  d)  =  [I  -f  P<f>(v,  -)]d. 

Thus,  equation  (1.3)  provides  a  complete  parameterization  of  the  class  of  all 
responses  {Ec}  that  can  be  obtained  by  internally  stable  control  of  the  system  E, 
with  the  stable  and  causal  system  (f>  serving  as  the  sole  parameter.  In  other  words, 
for  every  </>,  there  is  an  equivalent  controller  C  that  internally  stabilizes  £  and  yields 
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the  response  (1.3).  Conversely,  every  equivalent  controller  C  that  internally  stabilizes 
E  generates  an  explicit  (j>. 

On  a  superficial  inspection  of  (1.3),  the  attainable  performance  is  prescribed 
by  the  "numerator"  system  P  of  E  because  it  is  the  only  fixed  quantity  apart  from 
the  identity.  For  a  linear  time-invariant  system,  it  was  shown  in  [50]  that  the  at- 
tainable performance  was  determined  by  the  location  of  its  right  half-plane  zeros. 
For  a  nonlinear  system  E,  the  instabilities  of  P_1  (the  inverse  of  the  "numerator" 
system  P)  would  be  analogous  to  the  right  half-plane  zeros  of  a  linear  time-invariant 
system.  Using  this  analogy  it  could  be  inferred  that  the  instabilities  of  P"1  would 
limit  performance.  The  performance  bound  developed  in  the  dissertation  will  further 
confirm  this  hypothesis. 


1.2     Notation 
The  following  notational  convention  will  apply  unless  otherwise  stated. 

3£m  :  The  set  of  m-dimensional  real  vectors. 

S(^km)  :  The  set  of  m-dimensional  real  vector  sequences. 

Dx  :  The  A-step  shift  operator. 

£(E)  >  A  :  The  system  E  has  latency  of  at  least  A. 

|w|  :  The  £°°-norm  of  the  sequence  u. 

p(u)  :  The  weighted  f^-norm  of  the  sequence  u. 

||S||  :  The  Lipschitz  semi-norm  of  the  system  E. 

B  :  The  set  of  all  stable,  causal  and  recursive  systems  S:  T>  — >  S($lp] 

where  V  is  a  bounded  subset  of  S{$tm). 
Lip{V,S{W))  :  The  subset  of  8  with  ||E||  <  oo. 

llSy^jp  :  The  Lipschitz  norm  of  the  system  E. 

QmAu->^2)  '■  The  distance  from  any  u  6  <S"i  and  the  set  SV 

ge{Si,S-2)  '■  The  maximal  distance  from  any  u  G  Si  and  the  set  £2- 

V(Yj)  :  The  right  singularity  of  measure  of  the  system  E. 


CHAPTER  2 
TERMINOLOGY  AND  BASICS 

This  chapter  contains  a  summary  of  the  principal  mathematical  results  which 
are  needed  in  the  dissertation.  The  presentation  is  for  discrete-time  time-invariant 
nonlinear  systems.  Almost  all  the  necessary  mathematics  are  contained  in  [9,  14,  17, 
18,  22,  27]. 

The  set  of  real  numbers  is  denoted  by  3ft.  The  set  of  m-dimensional  real 
vectors  is  denoted  by  3ftm.  The  set  of  all  sequences  u  is  denoted  by  Sffi™)  where 
u  =  {uQ,  Ui,  ix2, . . .}  of  m-dimensional  real  vectors  ut  €  3£m,  i  =  0,1,2,....  Given 
the  sequence  u  6  S($tm),  the  ith  element  is  denoted  by  U{.  The  set  of  elements 
{«;,  i/j+i, . . . ,  Uj}  where  j  >  i  >  0  is  denoted  by  uj.  A  system  E,  from  an  input /output 
perspective,  is  a  map  E:5(3?m)  — »  5(5?p)  which  transforms  input  sequences  of  Tri- 
dimensional real  vectors  into  output  sequences  of  p-dimensional  real  vectors.  The 
image  of  a  subset  S  C  S(^R.m)  through  the  system  E  is  denoted  by  S[5].  The  entire 
image  of  system  the  £  is  denoted  by  ImE  where  Im  E  =  £[5(3ftm)].  Given  a  system 
£:5(3ftm)  — >  5(3^)  and  an  input  sequence  u  6  5,(3?m),  we  denote  by  Eu];  =  yt  the 
?'th  element  of  the  output  sequence  y  =  E«,  and  by  y  =  Sm]|  the  set  of  elements 
{yt,yl+u. . . ,  y3)  where  j  >  i  >  0  are  integers. 

There  are  two  kinds  of  binary  operations,  addition  and  composition.  For  a  pair 
of  systems  E1,E2:5(3?m)  — »  S'(3?p),  the  swm  is  defined,  as  usual,  by  (£i  +  S2)w  = 
Eiu  +  S2m  for  all  sequences  u  €  5(3ftm);  the  right  side  of  the  last  formula  is  the 
usual  elementwise  addition  of  sequences  of  real  vectors.  Composition  is  the  usual 
composition  of  maps. 
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Definition  2.0.1  A  system  S:S'(3?m)  — »  5(3ftp)  is  causal  (respectively,  strictly  caus- 
al,) if  it  satisfies  the  following  condition.  For  every  integer  i  >  0  and  for  every  pair 
of  input  sequences  u,v  £  5(3£m)  satisfying  ul0  =  vl0,  the  output  sequences  satisfy 
£w]o  =  £i>]o  (respectively,  Su]q+1  =  Tjv]1^1  ). 

A  system  M:5(3£m)  — >  S($lm)  is  bicausal  if  it  is  causal  and  if  it  possesses  a  causal 
inverse. 

A  system  E:  5(3£m)  — >  5'(9£p)  is  called  a  recursive  system  if  there  is  a  pair  of 
integers  rf,fji>  0  and  a  function  /:  (S^)^1  x  (3£m)^+1  ->  &p  such  that,  for  every  input 
sequence  t*  £  5(3£m),  the  corresponding  output  sequence  ?/  =  Ew  can  be  computed 
recursively  in  the  form 

Vk+v+i  =  f{yk,--  .,yk+v,uk,. . .  ,Ufc+M)  (2.1) 

for  all  integers  A:  >  0.  The  initial  conditions  y0, . . .  ,yv  must  of  course,  be  specified 
and  fixed.  The  function  /  is  called  a  recursion  function  of  E.  For  causality  (strict 
causality),  the  condition  r/  +  l  >//(//>  fi)  is  placed  on  system  E.  The  class  of  strictly 
causal  systems  include  every  system  E:5(5?m)  — >  S($lp)  that  can  be  represented  in 
the  form 

xk+1     =    f{xkluk) 

yk     =    h{xk),  k  =  1,2,....  (2.2) 

Here,  u  £  5(3£m)  is  the  input  sequence;  y  £  S'(3?p)  is  the  output  sequence;  and 
.r  £  S($ln)  is  an  intermediate  sequence  of  "states."  In  the  case  that  the  maps  /:  3ftn  x 
3£m  — >  3£n  and  /?.:  3£n  — >  3£p  are  continuous,  then  the  system  (2.2)  constitutes  a 
continuous  realization  of  the  system  E. 

For  a  real  number  9  >  0,  the  set  of  all  vectors  in  3£m  with  components  in  the 
closed  interval  [— 0,0]  is  denoted  by  [— 919)m.  The  set  of  all  sequences  u  £  5(3£m) 
with  elements  uz  belonging  to  [— 9,9]m  for  all  integers  i  >  0  is  denoted  by  S(9m). 
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Thus,  S(9m)  consists  of  all  sequences  bounded  by  0.  It  follows  then  that  a  system 
E:  S($m)  ->  S{W)  is  BIBO  (Bounded-Input  Bounded-Output)-s*afe/e  if  for  every  real 
number  0  >  0,  there  exists  a  real  number  M  >  0  such  that  £[,S(6>m)]  C  S{MP).  A 
sequence  u  E  S(3ftm)  is  said  to  be  bounded  if  there  is  a  real  number  6  >  0  such  that 

u  e  S(6m). 

The  basic  notion  of  stability  that  is  used  in  this  dissertation  is  related  to  con- 
tinuity with  respect  to  a  metric.  Two  norms  are  particularly  useful  in  this  context  to 
derive  a  metric:  the  f°°-norm  and  the  weighted  ^°°-norm.  The  f°°-norm  is  denoted  by 
|  •  |;  for  a  vector  a  =  (al5  a2, . . .  ,am)  E  3£m,  it  is  simply  \a\  =  max{|ai|,  |a2|, . . . ,  |am|}. 
For  a  sequence  u  E  5(3£m),  the  ^°°-norm  is  given  by 

\u\  =  sup  \ui\.  (2.3) 

The  weighted  ^°°-norm  is  denoted  by  pt  and  is  given  by 

PeM^supfl  +  eJ-'W,     e>0  (2.4) 

i>0 

for  a  sequence  u  E  5(3?m).  For  purposes  of  this  dissertation  different  values  of  e  will 
not  affect  our  results;  hence,  the  subscript  will  be  dropped  and  the  weighted  f°-norm 
is  simply  denoted  by  p.  The  use  of  the  weighted  f^-norm  simplifies  mathematical 
arguments  over  the  f°°-norm  because  the  bounded  set  of  sequences  S(6m)  is  compact 
with  respect  to  p.  The  norm  p  induces  a  metric  p  on  a  given  S'(3?m),  for  every  pair 
of  elements  w,  v  E  5'(§?m),  by  p(u,v)  =  p(u  —  v).  Formally,  the  notion  of  stability 
employed  in  this  dissertation  is  as  follows. 

Definition  2.0.2  A  system  E:S'(3?m)  — >  S'(3?p)  is  stable  with  respect  to  the  metric 
p  if  it  is  BIBO-stable,  and  if  the  restriction  S:5(am)  — >  S{$tp)  is  continuous  with 
respect  to  p  for  every  real  number  a  >  0. 

Definition  2.0.2  is  usually  referred  to  as  input/output  stability.  The  following 
concept,  which  describes  a  weak  form  of  uniform  continuity  with  respect  to  the  i°°- 
norm,  plays  a  fundamental  role  in  stabilization  theory  (see  [17]). 
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Definition  2.0.3  A  stable  system  Z:S($m)  ->  S(W)  is  differentially  bounded  if 
there  is  a  pair  of  real  numbers  e,  9  >  0  such  that,  for  every  pair  of  sequences  u  G 
S(&m)  and  v  G  S{em),  one  has  |E(u  +  v)  -  S(u)|  <  0. 

So  far,  only  stability  properties  of  individual  systems  have  been  mentioned. 
When  several  individual  systems  are  combined  into  a  composite  system,  a  stronger 
notion  of  stability  is  required,  and  it  is  usually  referred  to  as  internal  stability.  Internal 
stability  guarantees  desirable  stability  properties  of  the  composition,  and  takes  into 
account  the  effects  of  various  disturbances  and  noises  that  may  affect  the  component 
systems.  Consider  a  composite  system  E^  that  consists  of  s  individual  systems, 
labeled  S1, . . . ,  Ss  where  E*':  S(^m^)  ->  .S'(^p(j)),  i  =  1,. . .  ,s.  Individual  entries  in 
the  list  S1,...,!!'  may  represent  summers,  multipliers,  etc.  Let  u  G  5(3?m)  be  the 
external  input  sequence  of  the  composite  system,  and  let  y  G  S(3lp)  be  its  output 
sequence.  Let  uJ  G  5'(5?m(-7))  be  the  input  sequence  of  the  system  EJ  within  the 
configuration,  and  let  yJ  G  5(5^-^)  be  its  output  sequence.  The  interconnections 
among  the  subsystems  are  then  characterized  by  a  set  of  equalities  ul  =  y^l\  which 
determine  to  which  output  each  input  is  connected.  The  external  signal  u  is  now 
augmented  by  s  new  input  signals  rf  G  Sffi171^),  i  =  1, ...  ,5,  and  set  ul  =  y3^  -\-nl. 
For  each  i,  the  rf  acts  as  an  additive  disturbance  on  the  input  port  of  the  system  E\ 
The  disturbances  are  all  assumed  to  be  bounded  by  a  real  number  6  >  0,  so  that  in 
fact  rf  G  S(<Sm(!)),  i  =  l,..., 5. 

Let  E*s*:S'(ftm)  x  5(3?m<1))  x  •••  x  S{^m{s))  ->  S{W)  x  S{$p{1))  x  •••  x 
5(3ftp^):  (it,  J/1,  •  •  •  ,VS)  H^  ^""("i7/1:  •  •  •  IVs)  denote  the  system  induced  by  the  in- 
terconnected system  E^  and  the  disturbances,  having  the  input  signals  u,  r/1, . . . ,  ns 
and  the  output  signals  y,  y1, . . . ,  ys,  respectively. 

Definition  2.0.4  The  composite  system  S^  is  internally  stable  if  the  system  E*s*  is 
stable  in  the  sense  of  Definition  2.0.2.  The  composite  system  E^  is  strictly  internally 
stable  if,  besides  being  stable,  the  system  E*s*  is  also  differentially  bounded. 
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Definition  2.0.5  A  system  Y,:S{?Rm)  — >  5(3£p)  is  entirely  stabilizable  if  there  is  a 
strictly  internally  stable  control  configuration  that  stabilizes  S  over  the  entire  input 
space  S{$m). 

2.1      Right  Fraction  Representations  and  Coprimeness 

A  right  fraction  representation  of  a  system  H:S(^R.m)  — >  S'(3?p)  is  determined 
by  three  quantities:  a  subset  S  C  S(ffl),  q  >  0,  called  the  factorization  space;  and 
two  stable  systems  P:  S  — ►  S($p)  and  Q:  S  — >  5(9£m),  where  Q  is  invertible,  such 
that  £  =  PQ'1 .  A  right  fraction  representation  S  =  PQ~l  is  coprime  whenever 
the  stable  systems  P  and  $  are  right  coprime  according  to  the  following  definition 
([16,  18]).  (Let  G:  Sx  ->  52  be  a  map,  where  St  C  S{$m)  and  S2  C  S(3fcn)  are  subsets. 
For  a  subset  5  C  5(3£n),  we  denote  G*[S]  the  inverse  image  of  S  through  G,  i.e.,  the 
set  of  all  sequences  «G  5i  satisfying  Gu  £  5.) 

Definition  2.1.1  Lei  5  C  S(3fc9)  be  a  subset.  Two  stable  systems  P:  S  ->  S(9£p)  and 
Q:  5"  — >  .S"(^m)  are  right  coprime  whenever  the  following  conditions  hold. 

1.  For  every  real  number  r  >  0  there  is  a  real  number  9  >  0  such  that 

P*[S{Tp)}nQ*[S{Tm)]  c  S(0q). 

2.  For  every  real  number  r  >  0  the  set  S  fl  S(rq)  is  a  closed  subset  of  5(r9)  fu;?Y/i 
respect  to  the  topology  induced  by  p). 

The  concept  of  a  homogeneous  system  is  of  key  importance  to  the  theory  of 
right  coprime  fraction  representations  of  nonlinear  systems.  A  homogeneous  system 
has  the  property  of  being  a  continuous  map  whenever  its  outputs  are  bounded.  The 
precise  definition  is  as  follows. 

Definition  2.1.2  A  system  E:5(5£m)  — »  S(^p)  is  a  homogeneous  system  if  the  fol- 
lowing holds  for  every  real  number  a  >  0:  for  every  subset  S  C  S(am)  for  which 
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there  exists  a  real  number  r  >  0  satisfying  £[5]  C  S(tp),  the  restriction  of  £  to  the 
closure  S  of  S  in  S(am)  is  a  continuous  map  £:  S  — ►  S(tp). 

The  importance  of  homogeneous  systems  to  the  dissertation  is  stated  in  the 
next  two  theorems. 

Theorem  2.1.3  [17]  An  infective  system  £:S(3ftm)  — ►  5(3£p)  has  a  right  coprime 
fraction  representation  if  and  only  if  it  is  a  homogeneous  system. 

Theorem  2.1.4  [17]  Let  £:  S(3£m)  -»  5(3RP)  be  a  recursive  system.  7/E  has  a  recur- 
sive representation  yk+v+i  =  f{yk,  •  •  • ,  yjfc+rj,  w^, . . . ,  Ufc+M)  with  a  continuous  recursion 
function  f ,  then  £  is  a  homogeneous  system. 

The  following  two  theorems  are  important  for  stabilization  theory.  Theo- 
rem 2.1.5  defines  the  technical  details  in  the  stabilizing  system  E  in  Figure  2.1  with 
7r  =  B~l  and  ip  =  A.  Theorem  2.1.6  uses  a  basic  property  or  right  coprime  factor- 
ization, which  is  the  denominator  system  contains  the  exact  information  about  the 
instabilities  of  the  system,  of  right  coprime  fraction  representation  for  stabilization. 

Theorem  2.1.5  [19]  Let  E:  S(am)  — >  S(^.p)  be  a  system  with  a  bounded  input  space 
S(am),  a  >  0,  and  assume  it  has  a  right  coprime  fraction  representation  E  =  PQ~l , 
where  P:  S  — >  S{$tp)  and  Q:  S  — ►  S(am),  and  where  S  C  5(3?9)  for  some  integer  q  > 


+  . 
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Figure  2.1:  The  block  diagram  of  E^ 
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0.  Then,  for  every  stable  system  M:  S  — >  S  with  a  stable  inverse  system  M~l\  S  — >  S, 
there  exists  a  pair  of  stable  systems  A:  S($p)  — ►  S(5lq)  and  B:S(am)  — »  S(^Stq)  such 
that  AP  +  BQ  =  M. 

Theorem  2.1.6  [22]  Let  E:5($Rm)  — >  5'(5RP)  6e  a  system  having  a  right  co-prime 
fraction  representation  E  =  PQ~l ,  where  P:  S  — ►  5(3ftp),  Q:  S  —>  S(dtm),  S  C 
5(3?m),  ana1  where  Q  is  bicausal.  Let  D:S($r)  — >  5(3£m)  6e  any  s^a6/e  ana'  causal 
system  for  which  the  composition  ED:  S(3£r)  — >  S(3£p)  is  stable.  Then,  there  is  a 
stable  and  causal  system  (f>:  (3Rr)  — *  S  such  that  D  =  Q<f>. 

In  Figure  1.1,  E:5(3ftm)  — »  S(3ftp)  is  a  strictly  causal  system  that  needs  to  be 
controlled.  It  will  be  convenient  to  regard  the  external  input  sequence  v  as  a  fixed 
"parameter"  while  regarding  the  disturbance  d  as  an  external  input,  i.e.,  to  consider 
appropriate  partial  functions.  Since  no  restriction  will  be  placed  on  v  and  a1,  this  will 
have  no  effect  on  the  validity  of  the  final  result.  In  the  same  spirit,  we  shall  use  the 
notation 

^{v)z  =  C(v,z) 

in  which  v  can  be  intuitively  viewed  as  a  parameter  of  the  system  \P(u),  while  z  is 
its  input.  Control  is  achieved  by  a  causal  nonlinear  dynamic  controller  C:S(3Rm)  x 
S{$p)  ->  S(&m):  (v,  z)  h->  C(t/,  z).  For  Figure  2.2,  V(v)z  =  C{v,  z)  was  defined  with 
the  system  $(u):  S(3fcp)  ->  S{$m):z  ^>  ty(v)z  =  u.    A  result  of  [22]  is  now  stated 


Figure  2.2:  The  block  diagram  of  E$(„). 
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formally  in  Theorem  2.1.7.  which  furnishes  a  parameterization  of  the  class  of  all 
systems  that  can  be  obtained  from  a  given  system  E  by  internally  stable  control. 

Theorem  2.1.7  [22]  Let  E:5(3ftm)  — >  S($lp)  be  a  strictly  causal  system  having 
a  right  coprime  fraction  representation  E  =  PQ~l  with  a  bicausal  denominator 
Q:S(^R.m)  — »  S($m).  Assume  that  S  can  be  strictly  internally  stabilized  by  a  con- 
troller that  admits  representation  Ci(s,z)  =  G~1(z)[s  +  Tz]  =  u.  Then,  referring  to 
(1.2),  the  following  is  true.  The  class  of  all  input&  disturbance /output  responses  Ec 
that  can  be  achieved  through  internally  stable  control  of  E  is  given  by 
{Zc(v,d)  =  [I  +  P(f>(v)}d,  </>(•):  S(&m)  x  S{W)  ->  S{$m):(v,d)  ^  <f>(v)d  is  a  stable 
and  causal  system}. 

A  special  case  of  Theorem  2.1.7  is  when  the  system  E  is  stable.  In  that  case, 
the  right  coprime  factorization  of  S  =  PQ~X  can  be  taken  as  P  =  E  and  Q  =  I ,  and 
the  following  is  obtained. 

Corollary  2.1.8  [22]  Let  S:S(9£m)  -»  S($lp)  be  a  strictly  causal,  stable  and  differ- 
entially bounded  system.  Then,  referring  to  (1.2),  the  following  is  true.  The  class 
of  all  input& disturbance/ output  responses  Ec  that  can  be  achieved  through  internally 
stable  control  of  E  is  given  by 

{Zc{v,d)  =  [I  +  E<f>(v)]d,  <f>{-):  S{$m)  x  S{$p)  ->  S(9lm):(v,d)  k-»  <f>{v)d  is  a  stable 
and  causal  system} . 

2.2     Generalized  Right  Inverse 

For  the  nonlinear  recursive  system  E:5'(3?m)  — >  S($lp)  to  possess  a  right  in- 
verse, the  system  E  is  required  to  be  surjective.  By  restricting  the  range  of  E  to  the 
image  of  E,  a  surjective  system  Er:  S(^Rm)  — *  ImE  is  obtained  with  its  right  inverse 
system  S*:ImE  -»  S($m).  Let  Z9:S(Rep)  -»  S{$m)  be  any  extension  of  E*  from 
ImE  to  the  whole  space  S(^tp).  Then,  for  every  element  j/  G  ImE,  it  is  evident  that 
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SE5?/  =  EE*y  =  y.  The  system  Es  is  a  generalized  right  inverse  of  the  system  E.  The 
generalized  right  inverse  of  system  E  is  non-unique  when  system  S  is  not  bijective. 
The  next  theorem  and  its  corollaries  from  [14]  state  that  a  recursive  system  has  a 
recursive  generalized  right  inverse. 

Theorem  2.2.1  [14]  A  recursive  system  T,:S{$tm)  — >  S{$lp)  has  a  recursive  gener- 
alized right  inverse  Z9:S{$P)  ->  S{$m). 

A  couple  of  notes  from  the  proof  of  Theorem  2.2.1.  Let  the  system  E*:  ImE  — »  5(3?m) 
be  the  recursive  system  represented  by 

uk+v  =  g(uk,...,Uk+n-i,yk,  •••,yfc+H-*H-i)-  (2-5) 

In  order  to  extend  the  domain  of  E*  from  ImE  to  all  S{W),  let  #e:  (9fcm)"  x  (%?)*>+ "+2 
— >  3ft771  be  any  extension  of  the  function  g.  Then,  the  recursive  system  E5:  S(Rep)  — > 
5(3?m)  represented  by 

Ufc+A*   =  ^e(Mfc,---,Wjt+^-l,^,---,?/fc+M+r;+l)  (2-6) 

is  a  generalized  right  inverse  of  E. 

Corollary  2.2.2  [14]  Let  E:5(5Rm)  — >  5(3ftp)  6e  a  recursive  isomorphism.  Then,  its 
inverse  E_1:S'(3?P)  — >  ^(S?"1)  zs  a/so  a  recursive  isomorphism. 

Theorem  2.2.1  cannot  be  generalized  to  the  case  of  an  arbitrary  domain  P. 
Nevertheless,  Theorem  2.2.1  can  be  generalized  to  the  case  when  the  domain  T>  is 
in  the  following  particular  form.  A  subset  T>  C  S(^,m)  is  recursive  if  there  exists 
an  integer  £  and  a  function  a  assigning  each  point  (ao, . . . ,  a^)  £  (3£m)^+1  a  subset 
<r(ao)  C  $lm  such  that  V  consists  of  exactly  all  sequences  u  €  S{$tm)  satisfying 
Uk+£+i  £  a(uk  )  for  all  integers  k.  The  function  a  is  called  the  generating  function 
of  Z>.  For  instance,  if  11:5(^9)  — »  5(3?m)  is  a  recursive  system,  then  Imll  is  a 
recursive  subset  of  S(^Rm). 
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Corollary  2.2.3  [14]     A  recursive  system  S:P  — >  S(^tp),   where  T>  is  a  recursive 
subset  of  5(3£m),  has  a  recursive  right  inverse  E*:ImS  —*  V. 


CHAPTER  3 
THE  PERFORMANCE  BOUND 


This  chapter  starts  by  discussing  causality  of  systems  and  Lipschitz  norms 
which  will  be  useful  in  establishing  the  existence  theorem  of  a  best  approximate 
right  inverse.  The  performance  bound  will  later  be  defined  as  a  lower  bound  for  the 
performance  index  of  the  approximate  right  inverse  optimization  problem. 

3.1     Causality  of  Systems 

Recall  a  causal  system  E  is  one  for  which  the  values  of  the  output  sequence 
£w  up  to  and  including  index  i  depend  only  on  the  values  of  the  input  sequence  u  up 
to  index  i.  This  was  precisely  defined  in  Definition  2.0.1. 

For  an  integer  A,  the  A-step  shift  operator  is  denoted  by  Dx ,  defined,  on  any 
sequence  u  by 

DH  =  «*-a  (3-1) 

for  all  integer  k  for  which  u^-x  exists. 

Definition  3.1.1  Let  T,:S\  — ►  S2  be  a  system,  where  S\  C  S(^R.m)  and  S2  C  S(^.p). 
The  system  E  has  latency  of  at  least  A  if  there  is  an  integer  A  such  that,  for  every 
pair  of  input  sequences  u,v  £  S\  and  for  every  integer  k  >  0,  the  equality  u0  =  v0 
implies  Su]J+A  =  Ev]J+A. 

We  write  £(E)  >  A  if  the  system  £  has  latency  of  at  least  A.  Intuitively, 
the  latency  represents  a  'time  delay'  incurred  in  the  propagation  of  changes  from  the 
input  of  E  to  the  output  of  E.  It  is  a  simple  consequence  of  the  definition  that  a 
system  E  is  causal  if  and  only  if  £(E)  >  0,  and  it  is  strictly  causal  if  and  only  if 
£(E)  >  1.  From  [17],  a  few  simple  properties  of  latency  are  listed. 
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Theorem  3.1.2  [17]  Let  E^Si  -»  S2  and  £2:S2  -*  S3,  where  S1  C  S{$m),  S2  C 
.DAl[S(3£p)]  ancf  S3  C  D1+2[S ($&)],  be  systems,  each  with  a  well-defined  latency, 
respectively  £(£1)  >  Ai  and  £(£2)  >  A2.  Under  this  hypothesis  the  composition 
£  =  £2£i  ^as  a  well-defined  latency,  and  £(£)  >  Ai  +  A2.  /n  particular  i/Ai+A2  >  0 
£/ien  £  is  a  causal  system. 

Theorem  3.1.3  [17]  Le*  S:  Si  ->  S2,  w/iere  Si  C  S{$m)  and  S2  C  S(^p),  be  a  re- 
cursive system  with  a  recursive  representation  yk+v+i  =  f{yki  ■  •  •  3  2/fc+rn  w^, . . . ,  Ufc+Ai) 
fanG?  some  fixed  initial  conditions).  Under  this  hypothesis  £  /ms  well-defined  latency, 
and  £(T,)  >  t]  +  I  —  /.i . 

Theorem  3.1.4  [17]  A  system  £:  Si  ->  S2,  w/iere  Si  C  S(^m)  aW  S2  C  S(9£p), 
has  well-defined  latency  if  and  only  if  there  exists  an  integer  A  such  that  the  system 
D     £:  S'i  — »  Z}_A[S2]  is  a  causal  system. 

The  system  £:  Si  — >  S2,  where  Si  C  S(3?m)  and  S2  C  S(3£p),  is  a  recursive 
system  with  a  recursive  representation  yk+v+i  =  f(yk,---,yk+v>uk,---,Uk+v.)-  By 
restricting  the  range  of  £  to  the  image  of  £,  we  obtain  the  map  £r:Si  — »  £[Si] 
which  is  evidently  surjective  and  possesses  a  right  inverse  £*:  £[Si]  — >  Si  such  that 
££*?/  =  y  for  y  £  £[Si].  From  Theorem  2.2.1,  the  system  £*  is  recursive  with  a 
recursive  representation  given  by 

uk+n  =  g{uk,-  ■  ■  ,uk+fl-i,yk,-  ■  ■  ,yk+v+n+i)>  (2-5) 

Examining  (2.5),  n  +  1  shifts  are  required  on  the  output  of  a  right  inverse  system  £* 
to  guarantee  causality,  i.e.,  £(Z)('?+1'£*)  >  0. 

Lemma  3.1.5  Let  £:  Sx  ->  S2,  where  Si  C  S(&m)  and  S2  C  S(3£p),  be  a  recursive 
system  with  a  recursive  representation  yk+^+i  =  /(z/jt,  •  •  • ,  2/ *:+??,  w^, . . . ,  Ujt+M).  £e£ 
£*:  £[Si]  — *  Si  6e  a  n'g/ii  inverse  of  £  fry  restricting  the  range  of  £  £0  £/ie  image  of 
£  u^/i  a  recursive  representation  Uk+n  =  g{uk,.  ■  ■ ,  w^+^-ij  ?//t,  •  •  •  ,  2/fc+M+^+i)-  Under 
this  hypothesis  the  system  £  /ias  a  latency  of  at  most  7/  +  1,  i.e.,  £(Z)^+1^£*)  >  0. 
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3.2     Systems  with  the  Lipschitz  Norm 

In  the  previous  chapter  we  discussed  two  norms:  the  ^°°-norm  and  the  weight- 
ed f^-norm.  We  will  need  to  define  a  new  norm  for  a  set  of  systems  that  has 
the  convergence  property  which  is  that  every  Cauchy  sequence  in  a  compact  metric 
space  {X,p)  converges  to  some  point  of  X.  The  Lipschitz  norm  for  a  set  of  systems 
is  introduced  here  because  it  has  the  necessary  convergence  property.  The  Lipschitz 
norm  will  be  derived  from  a  semi-norm.  The  semi-norm  definition  from  [41]  is  given 
below. 

Definition  3.2.1   A  semi-norm  on  a  vector  space  X  is  a  real-valued  function  p  on 
X  for  all  x  and  y  in  X  and  all  scalars  a  such  that 

1.  p{x  +y)  <  p{x)  +p{y) 

2.  p(ax)  —  \a\p{x) 

3.  p(x)  t^O  ifx  t^O. 

In  this  section,  we  let  (5(3£m),  p)  and  (5(3?p),  p)  be  two  metric  spaces  over  the 
set  of  all  sequences  of  m-dimensional  and  p-dimensional  real  vectors.  Let  B  be  the 
set  of  all  stable,  causal,  and  recursive  systems  E:P  — >  S{$tp)  where  T>  is  a  bounded 
subset  of  S(^m).  Introduce  an  operator  ||  •  \\:B  — >  3ft+  defined  by 

||E||    =  sup .  (3.2) 

ui,u2GV,ui^u2  P\U\  —  U2) 

The  number  ||E||  has  the  following  properties. 

Lemma  3.2.2  //  E,  ^  G  B  and  ae$,  then 

1.  ||E||  =  0  if  and  only  if  E  is  a  constant  system  on  T>; 

2.  llaEII  =  \a\  ||S||;  and 
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3.  ||E  +  #||  <  ||S||  +  || tt||. 

The  number  ||E||  is  called  the  Lipschitz  semi-norm  ([7])  of  the  system  E  on  V. 
A  semi-norm  has  the  property  that  ||E||  =  0  does  not  necessarily  imply  that  S  =  0. 
In  fact,  it  can  be  easily  seen  that  ||E||  =  0  if  and  only  if  E  is  a  constant  system  (need 
not  be  zero)  that  maps  all  sequences  from  T>  to  the  same  sequence  in  S(§tp).  The 
Lipschitz  semi-norm  is  like  a  derivative  bound  or  gain  for  all  systems  in  B.  Let's 
define  Lip(V,S{'$tp))  as  the  subset  of  B  satisfying  ||E||  <  oo. 

It  is  clear  that  an  element  E  of  B  is  in  Lip(T>,  S($lp))  if  and  only  if  there  is  a 
number  L  >  0  such  that 

/?(S(ui)  -  S(u2))  <  Lp(ui  -  u2) 

for  all  ui,u2  £  T>.  Moreover,  ||E||  is  the  "least"  such  number  L.  It  is  also  evident 
that  a  system  with  the  Lipschitz  semi-norm  is  both  bounded  and  continuous  on  its 
domain.  The  Lipschitz  semi-norm  will  be  used  to  define  the  boundary  of  a  set  from 
which  an  approximate  right  inverse  can  be  chosen.  The  semi-norm  ||  •  ||  can  be  made 
into  a  norm  as  seen  in  the  following  theorem. 

Theorem  3.2.3  Let  u0  be  an  element  of  T>  and  let 

||S||Ltp    d^f    ,(E(uo))  +  ||E|| 

/w      uj                           P(E("i)  ~S(m2))  ,qqx 

=     p(S(wo))+  sup ,  {i.i) 

u1,u2eV,ui^u2  P{Ul—U2) 

then  the  number  ||E||^-p  defines  a  norm  for  all  E  £  Lip(T>,  S($lp)). 

||  E  || Lip  will  be  called  the  Lipschitz  norm  ([7])  of  the  system  E  defined  by 
u0  £  V.  A  convenient  choice  of  u0  is  of  course  u0  —  0  if  0  £  V,  where  note  that  E(0) 
is  not  zero  in  general.  To  prove  the  theorem,  it  amounts  to  showing  that  ||E||ljp  =  0 
implies  E  =  0,  the  zero  system.  This,  however,  is  an  immediate  consequence  of 
part  1.  of  Lemma  3.2.2. 
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Theorem  3.2.4   [7]    The  family  Lip(V,S(?Rp))  is  a  complete  metric  space  under  the 
Lipschitz  norm  ||  •  \\up- 

A  complete  metric  space  (Lip(V,  5(3£p)),  ||  •  \\up)  has  the  property  that  every 
Cauchy  sequence  of  systems  {£&}  converges. 

3.3     Approximate  Right  Inverse 

In  [50],  the  concept  of  an  approximate  inverse  was  formulated  for  linear  sys- 
tems. The  problem  of  disturbance  attenuation  was  addressed  with  a  configuration 
very  similar  to  Figure  1.1.  As  a  result  of  [50],  the  approximate  invertibility  of  the 
system  was  shown  to  be  a  necessary  and  sufficient  condition  for  disturbance  attenu- 
ation. The  optimization  problem  from  [50]  used  the  "H^-norm  ||  •  ||oo  and  a  weighted 
semi-norm  ||  •  \\w  with  the  property  ||  •  \\w  <  ||  ■  ||oo-  The  weighted  semi-norm  has 
a  weighing  filter  where  W  will  denote  a  stable  and  causal  system  of  unit  norm  such 
that  \\d\\w  =  \\WdWoo. 

The  definition  from  [50]  of  an  approximate  inverse  is  as  follows.  The  set  of  all 
stable,  causal,  and  linear  systems  is  denoted  by  Cs-  For  any  stable  and  strictly  causal 
system  P,  an  approximate  right  inverse  of  P  is  any  stable  and  causal  system  $  for 
which  ||7  —  P\P||w  <  ||/||w-  The  right  singularity  measure  of  P,  denoted  by  //(P),  is 

KP)=M  I'7  ~P*\\w.  (3.4) 

In  general,  f.i{P)  is  a  number  in  the  interval  0  <  /j{P)  <  1. 

For  nonlinear  systems  we  cannot  define  a  weighted  semi-norm  with  the  with 
the  property  ||  •  \\w  <  ||  ■  ||oo-  As  a  result,  we  cannot  use  (3.4)  to  derive  a  right 
singularity  measure  of  E  between  0  and  1.  Nevertheless  we  can  define  an  approximate 
inverse  optimization  problem.  Let's  first  review  some  properties  of  right  inverses  for 
nonlinear  systems.  In  Section  2.2  the  generalized  right  inverse  was  introduced  for  a 
class  of  recursive  systems  in  the  form  of  (2.1).  For  systems  S  that  are  not  bijective, 
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the  generalized  right  inverse  of  systems  E  are  non-unique.  For  the  systems  that  are 
bijective,  the  generalized  right  inverse  is  stable  but  not  necessarily  causal. 

By  assumption,  the  system  E  is  strictly  causal.  This  implies  £(E)  >  1.  The 
design  integer  A,  where  A  >  0,  is  used  to  characterize  the  class  of  approximate  right 
inverse  systems  that  takes  in  account  the  latency  of  E.  The  design  integer  A  is 
selected  before  the  optimization  process.  In  the  three  examples  of  Chapter  4,  the 
design  integer  was  selected  at  the  minimum  latency  of  E,  i.e.,  when  £(E)  >  Ai,  the 
design  integer  A  =  X\. 

The  approximate  right  inverse  optimization  problem  of  the  stable,  strictly 
causal  and  recursive  system  E:D  — ►  5(3£p)  with  D,  a  compact  subset  of  S(3£m),  is 
shown  in  Figure  3.1  with  e  =  [I(u)  —  D~xY,ty(u)]  where  e  £  S(W).  The  performance 
index  is 

inf     sup     p{e)  (3.5) 

where  A  is  the  set  of  all  stable  and  causal  systems  $:  S(ap)  — >  S($m).  The  system  $ 
is  a  stable  approximate  right  inverse  of  E  while  the  best  approximant  $*  is  not  guar- 
anteed to  be  stable.  To  guarantee  stability  of  the  best  approximant  $*,  the  systems 
E  and  #  are  chosen  from  a  class  of  systems  Lip(T>,  5(3£p))  and  Lip(S(ap),  5(3?m)), 
respectively.  The  set  A  has  to  be  compact  subset  of  Lip(S(ap),  5(3?m)).  The  next  two 
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Figure  3.1:  The  approximate  right  inverse  optimization  problem  of  E:    inf  sup  /9(e). 
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paragraphs  will  explain  why  we  need  the  set  A  to  be  compact  to  guarantee  stability 
of  the  best  approximant  $*. 

Approximation  theory  is  needed  to  solve  the  problem  of  finding  a  best  ap- 
proximant ^*.  In  solving  the  best  approximation  problem,  additional  hypotheses  are 
placed  to  guarantee  the  existence  and  uniqueness  of  the  best  approximation.  Our 
practical  problem  is  that  we  need  approximate  inverses  that  are  stable,  causal,  and 
recursive  but  for  now  we  will  take  approximate  inverses  that  are  only  stable  and 
causal.  We  shall  be  therefore  be  interested  in  a  basic  existence  Theorem  3.3.1,  which 
gives  sufficient  conditions  to  guarantee  existence  of  closest  point. 

Theorem  3.3.1   [6]    Let  K  denote  a  compact  set  in  a  metric  space.   To  each  point  p 
of  the  space  there  corresponds  a  point  in  K  of  minimum  distance  from  p. 

In  the  previous  section,  the  class  of  systems  with  the  Lipschitz  norm  was 
introduced  so  that  we  can  create  a  compact  set  of  systems  in  a  metric  space.  Using 
the  hypothesis  of  Theorem  3.3.1  as  a  starting  point,  we  can  then  prove  the  existence 
of  best  approximate  right  inverse. 

3.3.1     The  Existence  Theorem  of  Best  Approximate  Right  Inverse 

For  clarity,  the  following  sets  are  defined  for  this  subsection  only.  Let  a  >  0  be 
a  real  number.  Denote  by  Lip(S(ap),  5(3£m))  the  set  of  all  stable  and  casual  systems 
$:S{ap)  ->  S(3ftm)  satisfying  ||$||  <  oo.  Let  V  be  a  compact  subset  of  S(Um)  and 
let  Lip(T),  S(^p))  be  the  set  of  all  stable  and  causal  systems  S:  V  — >  5(3£p)  satisfying 
||E||  <  oo.  Let  G  >  0  be  a  real  number.  We  denote  by  Ag  the  compact  subset 
Lip{S{ap),S{Wn))  given  by 

Aq  =  {V\V  6  Lip{S{ap),S{$m))  and  ||tf||  <  G]  .  (3.6) 

Theorem  3.3.2   Let  E  G  Lip{T>,  5(9^))  be  a  strictly  causal  and  recursive  system  for 
a  given  real  number  G  >  0  and  for  a  design  integer  A  >  0.    There  is  a  stable  and 
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causal  system  $*  G  Ag  such  that 

sup    P(l(u)-D-xZV*(u))  <     sup    p  (l(u)  -  D-Xm{u))  (3.7) 

uGS(aP)  u€S(qP) 

/or  all  $  £  AG- 

Proof:The  right  hand  side  of  (3.7)  is  a  continuous  function  of  $.  If  we  take  ^  =  0, 
then  supu€S/Qp)  p(I(u)  —  D~AE\P(u))  =  a.  There  exists  an  infimum  of  the  right  hand 
side  of  (3.7)  since  it  is  continuous  with  respect  to  ^  and  it  is  bounded  between  0  and 
a.  Let 

r=    inf      sup    p  (l(u)  -  ZTA£#(u))  ,  (3.8) 

and  suppose  that 


'^AG  uES(aP) 


sup    p  (l(u)  -  D-Xmk(uj)  =  rfcj      fc  =  0,1,2,..., 

ueS(aP) 

with  rfc  monotonically  decreasing  to  r.  Then,  for  fc  sufficiently  large 


(3.9) 


sup    p  (l(u)  -  ZrAE$fc(u))  <  r  +  1  <  M. 


(3.10) 


Since  the  sequences  of  systems  {/  —  D~xT,^k}  are  bounded  and  monotonic,  it  follows 
that  the  sequence  of  systems  {/  —  Z)_AE#*.}  converges.  This  implies  that  for  every 
e  >  0  there  exists  an  integer  TV  such  that  l,n  >  N  implies 


sup    p(\l(u)-D-xmn{uj\  -  [l(ti)-2rAEtfj(u)l)  <£  (3.11) 

ueS(aP)      VL  J        L  J/ 

or 

sup    p  (lTAEtfn(u)  -  ZTA£ty(u))  <  e.  (3.12) 

The  Lip(S(ap),  5(5RP))  systems  set  is  a  subset  of  the  stable  and  causal  sys- 
tems r:,S(ap)  — ►  S(^lp)  satisfying  ||r||  <  oo.  The  composition  of  two  systems  of 
bounded  Lipschitz  norm  is  also  a  system  of  bounded  Lipschitz  norm;  hence,  for 
all  $  G  Lip(S(ap),S{$tm))  and  a  given  E  G  Lip{V,  S(3fcp)),  the  composite  system 
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£#  £  Lip(S(ap),  S{$p)).  Now  we  define  fk  G  Lip(S{oiP),  S{$p))  such  that  fk  =  £^A. 
and  rewriting  (3.12)  we  have 

\\D-xfn-D-xfm\\Lip<£.  (3.13) 

This  establishes  lim^— oo  D~xfk  =  D~xf*  where  the  system  D~xf*  is  the  limit  of  the 
sequence  of  systems  {D~xfk}.  To  translate  D~xf*  into  a  Lip(S(ap),  S(^tp))  system 
we  get  Dx{D~xf*)  =  f*]f.  The  system  f*]f  has  its  first  A  outputs  y0,yi,  •  •  •  ,2/A-i 
set  to  zero  because  they  were  not  used  in  the  convergence  of  the  sequence  {D~xfk}-  It 
follows  by  the  completeness  of  Lip{S{ap),  S{$p))  that  the  system  f*]f  G  Lip{S{ap), 
S($p)). 

The  sequence  of  systems  {D~xfk}  converges  to  the  limit  system  D~  /*.  There- 
fore, the  sequence  of  systems  {D~xY,xl>k}  converges  also  to  the  limit  system  D~x f* . 
In  other  words,  for  every  e  >  0  and  for  all  u  G  S(am),  there  exists  an  integer  N  such 
that  l,n>  N,  this  implies 

sup    p  (lTA£tfn(u)  -  ZTAE#f(w))  <  e.  (3.14) 

ueS(aP) 

From  the  properties  of  a  system  of  bounded  Lipschitz  norm  we  can  we  write 

sup    p(D-xmn{u)  -  D-Xml{u)) 

u£S(aP) 

<     \\D-xZ\\Ltp     sup    p(¥n(u)-ty(u)).  (3.15) 

uES(aP) 

We  will  assume  ||D_aS||ljp  ^  0,  which  implies  that  S  is  not  the  zero  system.  The 
sequence  of  systems  {^k}  is  a  continuous  mapping  from  a  compact  metric  space 
(S(ap),p)  into  a  compact  metric  space  {T*,p).  Now  the  sequence  of  systems  {^k}  is 
a  sequence  in  a  compact  metric  space  (Ag,  ||  '  \\liP)-  Therefore  some  subsequence  of 
{^k}  converges  to  a  system  in  Ag-  Thus  we  have  limj_oo  ^k3  =  $*  where  $*  G  Ag 
for  the  subsequence  of  systems  {\Pfc,}. 

Now  we  can  define  a  subsequence  of  {D~xfk}  by  {D~xfkj}  such  that 
D~x fk    =  D~xE^kj-  Every  subsequence  of  {D~xfk}  converges  to  D~xf*  so 
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\'imj^00D~xfkj  =  D~xf*;  therefore 


>-a< 


>-A 


lim  D-*XV  =  D-*fm. 


(3.16) 


The  system  $*  G  Ag  is  best  approximate  right  inverse  of  the  system  £  with  optimal 
sequence  of  systems  {$k,}  €  Ag  that  converges  to  $*.      ■ 


3.4     The  Measure  of  Right  Singularity 

As  has  been  previously  noted,  for  every  pair  of  elements  u,v  €  S(3£m),  the 
metric  pt  is  defined 

pe(u,v)  =  sup  (1  +  e)""'  \u{  —  v,-|,     e  >  0.  (3.17) 

t>0 

The  distance  from  any  u  E  5(5?m)  and  the  set  S2  C  5(3£m)  is  defined  by  the  number 

def 


v€52 


(3.18) 


We  can  now  define  the  maximal  distance  from  any  u  £  S\  C  S(3£m)  and  the  set  S2 
by  the  number 

def 


Qc(S\,S2)  =  sup  qm((uiS2)=  sup    inf  pt{u,v). 
ueSi  ue5i  U£S2 


(3.19) 


From  Theorem  3.3.2,  there  exists  a  best  approximate  right  inverse  system 
$*  €  Ag  for  a  design  integer  A  >  0  such  that  the  inequality 

sup    pt  (l{u)-  D-Xm'{u))  <     sup    pAl{u)~  ZrA£tf(u))  (3.2O) 
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Figure  3.2:  The  performance  bound  optimization  problem  of  S:    supu€5(aP)  pc{e] 
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holds  for  all  ^  £  AG  and  for  all  e  >  0.  Using  $*  for  $  with  e  =  [I(u)  -  £>_A£#*(u)] 
in  Figure  3.1,  Figure  3.2  represents  the  performance  bound  optimization  problem  of 
the  causal  and  recursive  system  £  G  Lip(T>,  S(^.p)),  where  the  performance  index  is 

sup     Pe(e)  (3.21) 

u£S(aP) 

where  e  6  S($tp).  Now  we  want  to  put  a  lower  bound  on  the  term  supu€S(aP)  /?e(e)  = 
suPu€S(aP)  Pt(I{u)  ~  D~xT,$*(u)).  Using  (3.18)  and  (3.19)  we  can  write 

sup     pt  (l(u)  -  D-Xm*{u))     =       sup     6Mt(u,D-xm*[S{ap)]) 

ueS(aP)  V  '  u€S(aP)  V  ' 

=    Qt(S(ar),D-xm*[S(ap)}) 

>    ft(%p),fl-AImE).  (3.22) 

Note  that  even  if  ImS  is  unbounded,  ge  (S(ap),  Z)_AIm£J  is  bounded  since  the  first 
operand  of  Qe{-,-),  S(ap),  is  bounded. 

Definition  3.4.1  For  any  system  £  £  Lip(T),  S(^RP))  that  is  a  strictly  causal  recur- 
sive of  the  system  (2.1),  the  right  singularity  measure,  denoted  by  "P(£),  for  a  real 
a  >  0  and  a  design  integer  A  >  0  is 

P(£)  =f  ge  (S{ap),  ZrAImE)  .  (3.23) 

The  performance  bound  is  finally  defined  as  the  right  singularity  measure  "P(£). 
It  provides  a  measure  of  the  ability  to  match  a  set  by  subset  of  the  image  of  a  system 
as  seen  in  (3.23).  The  performance  bound  came  from  the  requirement  to  provide  an 
estimate  of  the  minimal  effect  of  the  disturbance  signal  d  on  the  output  signal  z  of 
Figure  1.1.  By  assumption  the  system  E:  V  — >  S'(3?p)  is  stable  where  T>  is  a  compact 
subset  of  5(3£m).  In  that  case,  the  right  coprime  factorization  of  £  =  PQ~l  can  be 
taken  as  P  =  £  and  Q  =  I  where  I:T>  —>  T>,  the  identity  system,  and  P:V  — >  £[P]. 

When  the  system  £:  S(am)  — >  S(W)  is  unstable,  the  right  coprime  factoriza- 
tion of  £  =  PQ-1  is  taken  as  P:  S  ->  E[S(am)]  and  Q:  S  ->  S{am)  where  5  C  S(9fem) 
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and  Q  is  bicausal.  From  Theorem  2.1.5,  any  stable  and  bicausal  system  M:  S  — >  S 
with  a  stable  inverse  system  M_1:  5  — >  5  can  be  selected,  so  there  exists  a  pair  of  sta- 
ble systems  A:S{W)  -»  S($m)  and  B:S(am)  ->  S(3£m)  such  that  AP  +  BQ  =  M. 
The  stabilization  of  £  can  be  seen  in  Figure  2.1  with  7r  =  B~x  and  <£>  =  A.  The 
closed  loop  response  of  the  stabilization  is  PM~X .  Now  M  can  be  selected  as  the 
identity  system  /:  S  — >  S.  The  factorization  space  S  can  be  taken  as  S(am)  and  [21] 
describes  the  construction  of  the  stabilizing  controllers  that  yield  the  factorization 
space  of  S(am).  As  a  result  the  closed  loop  response  is  P:  S(am)  — ►  T,[S(am)]. 

The  preceding  paragraphs  illustrate  that  the  "numerator"  system  P  directly 
drives  the  performance  bound  by  the  image  of  its  system  S.  The  image  of  E  for 
discrete-time  nonlinear  systems  is  analogous  to  the  right  half-plane  zeros  of  linear 
continuous-time  systems.  Both  of  these  quantities  are  fixed  and  cannot  be  altered  by 
feedback.  Now  we  can  use  the  performance  bound  for  systems  of  practical  origin. 


CHAPTER  4 
CALCULATION  OF  THE  PERFORMANCE  BOUND 


4.1     The  Estimate  of  the  Performance  Bound 

It  turns  out  that  the  calculation  of  the  right  singularity  measure  'P(E)  for 
nonlinear  systems  is  rather  laborious.  We  develop  in  this  section  an  estimate  "P(S) 
of  V(E)  which  is  easier  to  calculate.  This  estimate  satisfies 


P(E)<P(E). 


(4.1) 


From  Theorem  3.3.2,  there  exists  a  best  approximate  right  inverse  system  \P*  €  Ag 
such  that  the  inequality 

sup    p((l{u)-D-xW*{u))  <    sup    Pt(l{u)-D-xW{u))  (3.20) 

ueS(aP)        V  '         ueS(aP)        v  ' 

holds  for  all  \P  E  Ag  and  for  a  design  integer  A  >  0.  Theorem  3.3.2  demonstrates 
the  existence  of  a  best  approximate  right  inverse  system  but  does  not  provide  the 
construction  of  a  best  approximate  right  inverse  system. 


h, 

I 

*- 

w 

V 

u 

( 

v 

$ 

y 

4—1 

D~x 

1 

u  €  S(a 

n 

w 

w 

W 

e  i 

+ 


Figure  4.1:    The  estimate  of  the  performance  bound  optimization  problem  of  the 
system  E  with  an  approximate  right  inverse  system  $:    "P(E)  =  supuG5(QP)  /?f(e). 
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Using  a  causal  approximate  right  inverse  system  $  £  Lip(S(ap),  S(^m))  for  a 
given  design  integer  A  >  0  from  Figure  3.1,  Figure  4.1,  with  e  =  [I{u)  —  D_AE^(u)], 
represents  the  estimate  of  the  performance  bound  optimization  problem  of  the  causal 
and  recursive  system  £  £  Lip(D,S($lp)),  where  the  performance  index  is 

sup     p€(e)  (4.2) 

u£S{aP) 

where  e  £  5(3£p).  Now  we  want  to  express  the  term  supuGS(aP)  pt(e)  —  supueS,aP^  pf 
{I{u)  —  D~xT,^(u))  using  the  operator  ge{-,-)-  Using  (3.18)  and  (3.19)  we  can  write 

sup     Pi(l(u)-D-Xm(u))     =       sup     gMiUD-x^[S(ap)]) 

ueS(aP)  V  '  uES(aP) 


Qi(S(ap),D-xm[S(ap)}).  (4.3) 


Definition  4.1.1  For  any  system  £  £  Lip(V,  S($p))  that  is  strictly  causal  and  re- 
cursive, an  estimate  of  the  performance  bound,  denoted  by  "P(E),  given  real  a  >  0 
and  design  integer  A  >  0,  is  given  by 

-P(S)  d^f  qc  (S(ap),  D-xm[S{ap)])  ,  (4.4) 

where  $  £  Lip(S(ap),  S(9£m))  is  a  causal  approximate  right  inverse  system. 

4.2     Practical  Implementation  Techniques 

For  the  applications  presented  here,  we  will  be  interested  in  the  response  of 
the  systems  over  a  finite  interval  of  time,  say  for  a  discrete  set  [0,T],  where  T  >  0 
is  a  fixed  integer.  Let  Ct{^u)  be  the  set  of  all  functions  h:  [0,  T]  — »  3£n.  Denote  by 
\h\  =  supie[0iT]|/^)|,ther°-normon  CT{$n)  and  by  pe(h )  d=  supiG[0)T](l  +  e)-l\h{i)\, 
for  e  >  0,  the  norm  pt  on  CV(3£n).  Let  h  £  Ct(^u)  be  any  function.  Then 


and 


sup  (1  +e)~i\h(i)\  <   sup  \h{i)\ 

t'€[0,T]  t'€[0,T] 


(l+e)T  sup  (1  + e)-"|fc(t)|  =   sup(l  +  e)T->(i)|>    sup   \h(i)\, 

ie[0,T]  t'G[0,T]  *e[o,T] 


32 


so  that 


pt(h)  <  \h\  <  {1  +  e)1  pt{h)     and     (1  +  e)~'  \h\  <  pt(h)  <  \h\ 


(4.5) 


on  Cr{^n)-  These  inequalities  imply  the  equivalence  of  the  norms  |  •  |  and  pt(-)  on 
Cr(3ftn).  We  will  use  the  ^°°-norm  for  its  simplicity  in  the  calculation  of  the  estimate 
of  the  performance  bound  "P(E)  (so  the  norms  |  •  |  and  pt{-)  are  interchangeable  over 
finite  intervals  of  time). 

For  the  mathematical  analysis  presented  in  Chapter  2  and  Chapter  3,  the 
signals  of  interest  are  always  bounded.  For  implementation  purposes,  the  maxi- 
mum peak  deviation  of  the  input  signal,  max,>o  |Au,|  where  Aut-  =  u,+i  —  w,  is 
limited.  To  limit  the  maximum  peak  deviation,  the  signals  are  transformed  by  using 
a  discrete-time  lowpass  filter.  A  new  family  of  sequences,  denoted  by  SnLP{a),  is 
formed  by  using  a  second  order  discrete-time  lowpass  filter  with  a  cutoff  frequency 
Q,lp  in  rad/sec  on  S(a).  The  transformation  between  S(a)  and  SnLP(a)  is  bijective 
since  the  lowpass  filter  is  one-to-one  and  the  image  of  the  lowpass  filter  is  SnLP(a). 
Thus,  the  transformed  space  SnLP(a)  is  isomorphic  to  the  input  space  S(a). 

The  alternate  definitions  of  'P(Il)  are  given  in  (4.3)  and  using  the  £°°-norm 
(for  e  =  0)  we  get 


-P(S)  =     sup 

u£S(aP) 


I(u)-D-xXV(u)  =    sup     gMo(u,D-xm[S(ap)}).         (4.6) 

u&S{aP) 
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Figure  4.2:  The  Function  f:S{ap)  ->  3£+:  u  t-*  f(u). 
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Define  a  function  /:  S(ap)  — >  3?+  that  incorporates  the  lowpass  filter,  as  displayed  in 
Figure  4.2,  such  that 

£(£)=    sup     /(«),  (4.7) 

uGS(aP) 

where  the  function  /  is  defined  as  follows 

f{u)  =  \l{LPF{u))  -  D~x^  (LPF(u))\ .  (4.8) 

The  parameters  of  the  function  /  are  the  system  E,  an  approximate  inverse  system 
$,  the  design  integer  A,  the  optimization  set  S(ap),  and  the  lowpass  filter  cutoff 
frequency  Qlp-  The  function  /  must  also  satisfy  f(u)  =  qm0{u,  D~xT*ty[SnLP(ap)]) 
where  qm0{'i  ■)  is  defined  in  (3.18).  It  then  follows  that  /(•)  is  a  convex  function  when 
its  input  or  optimization  set  S(ap)  and  the  image  set  Y,ty[SnLP(ap)]  are  convex  since 
qm0  is  a  distance  function  from  a  point  in  the  optimization  set  to  the  image  set.  The 
distance  function  qm0  by  its  construction  is  a  convex  function.  In  practice,  the  image 
set  is  not  convex;  however,  studying  the  optimization  problem  of  finding  a  global 
maximum  of  a  convex  function  will  give  us  insight  to  the  non-convex  optimization 
problem. 

Unlike  the  minimization  of  a  convex  function,  once  a  local  maximum  has 
been  found  there  is,  more  or  less  by  definition,  no  local  information  to  tell  you  how 
to  proceed  to  a  higher  local  maximum.  In  particular,  there  is  no  local  criterion 
for  deciding  whether  a  given  local  maximum  is  really  the  global  maximum.  Hence, 
maximizing  a  convex  function  is  usually  a  much  harder  task  than  minimizing  a  convex 
function.  Theorem  4.2.1  shows  that  a  convex  function  achieves  a  maximum  over  a 
compact  polyhedral  set  at  an  extreme  point. 

Theorem  4.2.1  [3]  Let  g:  S{W)  — >  5ft  be  a  convex  function,  and  let  S  be  a  nonempty 
compact  polyhedral  set  in  S($lp).  For  the  problem  of  maximizing  g(u)  subject  to  u  E  S , 
there  is  an  optimal  solution,  u,  where  u  is  an  extreme  point  of  S  such  that  g(u)  >  g(u) 
for  all  u  G  S . 
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Using  Theorem  4.2.1,  he  global  maximum  of  f(u)  maybe  found  by  evaluating 
all  the  extreme  points  of  S(ap)  and  finding  the  maximum  on  this  set.  For  example, 
if  we  have  a  one-dimensional  set  of  sequences  which  is  a  finite  length  of  250  elements, 
the  number  of  extreme  points,  and  hence  the  number  of  evaluations  of  /,  would  be 
2250  =  1.8903  x  1075.  This  procedure  is  laborious  as  far  as  computation  is  concerned. 
Although  the  image  set  S  is  non-convex,  the  property  that  the  optimization  set  is 
convex  is  enough  to  suspect  that  global  maximum  would  be  found  on  an  extreme 
point  of  the  optimization  set.  If  we  take  an  interior  point  Ui  from  the  optimization 
set  which  is  a  nonempty  compact  polyhedral  set  and  form  a  neighborhood  NTi(ui) 
of  radius  r?,  then  we  can  find  a  m!+1  £  NTl(ui)  such  that  QM0(ui+i,  S)  >  QM0{un  S) 
and  uz+i  is  closer  to  the  boundary  than  ut.  The  radius  rt  of  a  neighborhood  Nri(ut) 
is  selected  so  that  we  do  not  accumulate  in  an  interior  local  maximum.  Since  the 
optimization  set  is  convex  and  qm0  is  a  distance  function,  the  sequence  {ut}  will 
converge  to  a  point  on  the  boundary.  Therefore,  the  global  maximum  will  be  on 
the  boundary.  For  simulation  purposes,  the  global  maximum  is  assumed  to  be  an 
extreme  point. 

The  three  systems  of  practical  origin  that  have  been  selected  for  investigation 
are  all  fourth  order  nonlinear  systems.  Each  of  the  systems  has  a  distinct  memoryless 
nonlinear  subsystem.  These  subsystems  represent,  respectively,  dead-zones,  position 
constraints,  and  rate  constraints  for  the  problems:  permanent  magnet  stepper  motor 
model  example,  multivariable  process  control  model  example,  and  aerodynamic  model 
example.  Nonlinear  programming  is  used  to  calculate  the  estimate  of  the  performance 
bound.  The  algorithms  that  are  developed  take  advantage  of  the  traits  of  each  of 
the  nonlinear  systems  and  of  the  finite  bandwidth  of  the  input  signal  space.  A  useful 
trait  of  a  system  is  a  property  that  would  help  reduce  the  region  of  optimization 
in  the  search  for  a  global  maximum.  It  will  be  shown  in  the  next  sections  that  the 
computation  is  not  as  laborious  as  previously  conjectured. 
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4.3     Permanent  Magnet  Stepper  Motor  Model  Example 

This  section  is  concerned  with  voltage  control  of  a  permanent  magnet  (PM) 
stepper  motor.  This  example  comes  from  [5],  wherein  the  authors  developed  a  model- 
based  control  law  using  exact  linearization  implemented  on  an  industrial  setup. 

The  PM  stepper  motor  model  consists  of  a  slotted  stator  with  two  phases 
and  a  permanent  magnet  rotor.  One  side  of  the  rotor  is  a  north  pole,  and  the  other 
side  is  a  south  pole.  The  teeth  on  each  side  of  the  rotor  are  out  of  alignment  by  a 
tooth-width.  The  equations  of  motion  for  the  PM  stepper  motor  model  are  given  by: 

~     =     j[va-Ria  +  KmUsm(Nr6)}  (4.9) 

^    =    j[vb-Rib-Kmucos(Nr0)]  (4.10) 

—    =    -  [-Kmia  sm(Nr6)  +  Kmib  cos(Nr0)  -  Bu]  (4.11) 


d6_ 
~dt 


=    u  (4.12) 


where  va,  vb  and  ia,  ib  are  the  voltages  and  currents  in  phases  A  and  B,  respectively. 
Further,  u  is  the  rotor  (angular)  speed,  6  is  the  rotor  (angular)  position,  B  is  the 
viscous  friction  coefficient,  J  is  the  combined  inertia  of  the  motor  and  the  translation 
table,  Km  is  the  motor  torque  constant,  R  is  the  resistance  of  the  phase  winding,  L  is 
the  inductance  of  the  phase  winding,  and  jVr  is  the  number  of  rotor  teeth.  Table  4.1 
shows  the  values  of  the  stepper  motor  parameters  used  in  [5]. 

Table  4.1:  The  PM  Stepper  Motor  Model  Parameters 


Parameter 

Value 

L 

1.5  mH 

R 

0.55  n 

J 

4.5  x  10"5  kg-m2 

Km 

0.19  N-m/A 

Nr 

50 

B 

8.0  x  10~4  N-m-s/rad 
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For  the  PM  stepper  motor  model,  there  is  an  appropriate  nonlinear  coordinate 
transformation  which  is  known  as  the  direct-quadrature  (DQ)  transformation.  The 
DQ  transformation  for  the  phase  voltages  and  currents  is  defined  as  follows: 


Vd 

def 

. v* . 

id 

def 

.  li  . 

cos(Nr6)  sin(JVr0) 

-sm(Nr$)  cos(Nr0) 

cos(Nr9)  sm{Nr0) 

-sm{Nr6)  cos(Nr0) 


Va 
Vb 

ia 


(4.13) 
(4.14) 


The  direct  current  id  corresponds  to  the  component  of  the  stator  magnetic  field 
along  the  axis  of  the  rotor  magnetic  field,  while  the  quadrature  current  iq  corresponds 
to  the  orthogonal  component.  The  application  of  the  DQ  transformation  to  the 
original  system  equations  (4.9)  through  (4.12)  yields: 


did 
dt 
diq 

dt 
du> 

~dt 
dfi_ 

dt 


—  [vd  -  Rid  +  NrLu>iq] 

—  [vq  -  Rig  -  NrLivid  -  Kmu] 


J 


[Kmiq  -  Bu] 


(4.15) 
(4.16) 
(4.17) 
(4.18) 


where  Vd  is  the  direct  voltage,  vq  is  the  quadrature  voltage,  id  is  the  direct  current, 
iq  is  the  quadrature  current,  u  is  the  angular  velocity,  and  6  is  the  angular  position. 
The  quadrature  component  iq  of  the  current  produces  torque  while  the  direct 
component  id  does  not  produce  any  torque.  However,  in  order  to  attain  higher  speeds, 
it  is  necessary  to  apply  a  negative  direct  current  in  order  to  cancel  the  effect  of  the 
back-emf  of  the  motor  which  results  from  the  "back-emf"  term  Kmu  in  (4.16).  The 
direct  current  id  will  be  forced  to  be  negative  by  the  correct  choice  direct  voltage  Vd- 
The  optimal  lead-angle  approach  is  used  in  [5]  to  properly  choose  Vd  which  becomes 
a  nonlinear  function  of  the  states.  The  equations  of  motion  are  discretized  using  first 
order  Euler  integration.    Let  xk  —  [idk  iqk  a?*  Ok]'.    Equations  (4.15)  through  (4.18) 
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take  the  following  form: 


Xk+i     =    f{xk)+g{xk)-v 


uk 


0    0    10 


<?Jk 


•  a:*, 


(4.19) 
(4.20) 


where  x^  €  5R4,  and  /,  g  are  vector  fields  on  3ft4. 

The  analytical  expression  for  the  input  dead-zone  is 

'  vCqk  -  1      if     vCqk  >  1 
vqk  =  DZ(vCqk)  =  <    0  if     -KvCqk<l 

k  vCqk  +  1      if     vc,fc  <  -1. 


(4.21) 


The  system  E:X>  — ►  5(3?)  :  vc9  l—>  w,  where  V  C  S(9£),  is  the  composition  of 
the  equations  of  motion  with  the  dead-zone  for  which  vcq  is  the  quadrature  voltage 
command.  Figure  4.3  shows  the  block  diagram  of  the  system  E. 

4.3.1     Simulation 

The  signal  of  interest  is  the  angular  velocity  signal  u.  It  is  assumed  that  all 
signals  of  interest  start  at  rest  (i.e.,  ujn  =  0  for  all  n  <  n0).  In  one  of  the  applications 
presented  in  [5],  the  desired  speed  trajectory  required  the  motor  to  go  from  to  zero 
to  1350  RPM  in  10  ms  with  a  sample  rate  of  10  kHz.  The  bandwidth  of  u  was  set  to 
50  Hz.  The  rise  time  was  about  10.6  ms  and  the  absolute  bound  of  \u\  <  1350  RPM 
was  imposed. 


Dead-Zone 


VCqk 


/ 

/ 

vq, 


\+,  =  f(x  J  +  g(x  „)  vq 


V      ^k 


ujk  =  [0  0  1  0]  x 


^jt+i 


V 


Figure  4.3:   The  block  diagram  of  the  system  E:X>  — >  5(3ft)  :   vcq  >->  uj  for  the  PM 
Stepper  Motor  Model  Example. 
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The  block  diagram  of  an  approximate  right  inverse  system  $:  S2n-so{a)  — ► 
5(3?)  :  ujc  l— ►  vqmyf  is  shown  in  Figure  4.4.  The  first  block  calculates  the  raw  quadra- 
ture voltage  v'qk  for  a  given  angular  velocity  command  uck-  The  analytical  expression 
for  the  continuous  dead-zone  inverse  (CDI)  is  given  by 


U,invk  =  CDI(v'  )  =  i 


r<  +  i    ^    <>o± 


G-v' 


q\ 


if  -g=t<u;*<g=t 


(4.22) 


<"1      if     «k<-A 


Here  G  =  \\CDI(-)\\  represents  the  Lipschitz  semi-norm  of  the  CDI  system.  The 
range  of  G  is  1  <  G  <  oo.  In  the  CDI  block,  the  voltage  quadrature  command 
vqinvk  is  determined  by  using  (4.22)  on  the  raw  quadrature  voltage  v'qk.  The  voltage 
quadrature  command  vqmVk  is  then  used  to  update  the  states  xmVk  of  the  system  $. 

Two  integrations  of  vq  are  required  to  reach  u>.  This  results  in  the  discretized 
model  of  at  least  a  two  step  delay.  Therefore,  £(E)  >  2.  From  Figure  4.4  it  follows 
that  £(ty)  >  0.  Theorem  3.1.2  applied  to  the  composite  system  E$:  52*. 50(a)  — ► 
5(3?)  :  uc  ■->  w  yields  £(£$)  >  2.  For  the  simulation,  g(xk)  =  [0  ±A*  0  0]'  ^  0  in 
(4.19)  for  all  vcq  £  T>  and  a  given  x0.  The  design  integer  is  taken  as  A  =  2  for  an 
approximate  right  inverse  system  $. 

For  the  simulation  of  40  milli-seconds,  a  sampling  period  of  100  micro-seconds 
yields  a  sequence  of  400  points.    In  the  previous  section,  the  global  maximum  was 


idk 


rxk=f(xk.1)  +  g(xk.1)vqk  ! 

I 
\  cJfc  =  [0  0  1  0]  x    k  I 


-vq 


Continuous 
Dead-Zone  Inverse 


Dead-Zone 


xk  =  f(xk.I)  +  g(xk.|)vqk     |*vq 


7" 


/ 


$ 


<7invfc 


Figure  4.4:  The  block  diagram  of  an  approximate  right  inverse  system  $:  52ti-5o(c*! 
5(3?)  :  uc  ^— ►  fqinv  for  the  PM  Stepper  Motor  Model  Example. 
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shown  to  be  on  the  boundary.  For  simulation  purposes,  the  global  maximum  is 
assumed  to  be  an  extreme  point.  The  pre- filtered  binary  random  sequences  have 
only  the  values  ±1350  RPM,  the  extreme  points  of  5(1350).  Table  4.2  shows  the 
Monte  Carlo  results  for  4  different  system  gains  G  of  the  GDI  system:  2,  10,  100 
and  1000,  where  the  error  \u>c  —  D~2Y,ty(u>c)\  is  evaluated  for  each  of  the  sequences. 
As  can  be  seen  from  Table  4.2,  as  G  — >  oo,  the  error  \uc  —  D~2Y,ty(uc)\  — >•  0  in 
a  linear  fashion.  For  the  limiting  case  as  G  — ►  oo,  an  approximate  right  inverse 
system  ^  is  not  continuous.  However,  the  performance  bound  'P(E)  appears  to  be 
directly  proportional  to  the  system  gain  G.  This  observation  will  be  proven  true  in 
the  following  subsection. 

From  Table  4.2,  the  maximum  value  of  \uc  —  Z)~2S\I>(u;c)|,  for  G  =  10  rad/sec, 
is  0.02686.  This  corresponds  to  the  18th  sequence  in  the  Monte  Carlo  run.  In  the 
following  subsection  we  will  calculate  the  estimate  of  the  performance  bound  V(T,) 
over  the  optimization  set  52*50(1350).  A  first  estimate  is  "P(S)  >  0.02686.  In 
Figure  4.5,  the  18th  sequence  for  wp,  the  filtered  square  wave,  is  plotted  with  the 
response  sequence  Ev^^c),  without  the  D~2  shift  since  the  errors  are  less  than  0.01% 
of  the  magnitude  of  the  signal.  The  absolute  angular  velocity  error  sequence  is  plotted 
in  Figure  4.6.  The  quadrature  voltage  vq  and  the  absolute  quadrature  voltage  error, 
\v'qk  —  vqk\,  are  plotted  in  Figure  4.7  and  Figure  4.8,  respectively.  Note,  that  Figure  4.6 
and  Figure  4.8  appear  to  be  the  same  figure  but  of  different  scales. 

4.3.2     The  Estimate  of  the  Performance  Bound  V(T,) 

In  this  section,  we  will  find  an  estimate  of  the  performance  bound  V(E)  for  the 
voltage  control  problem  of  a  PM  stepper  motor  by  using  the  optimization  definition 
as  follows: 

£(£)=  sup  uc  -  D~2m(ujc)  (4.23) 

wc€S2*-.5o(1350) 
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Table  4.2:  Family  of  Sequences, 

S27r.5o(1350),  Comparison 

Sequence 

Approximate  Right  Inverse  Error  \ojc  ~  j 

D-2Etf(tjc.)l 

Gain  G  —  2 

Gain  G  =  10 

Gain  G  =  100 

Gain  G  =   1000 

1 

0.1339 

0.02628 

0.001151 

0.0002244 

2 

0.1341 

0.02588 

0.001021 

0.0000000 

3 

0.1342 

0.02510 

0.001979 

0.0001157 

4 

0.1343 

0.02546 

0.001346 

0.0000000 

5 

0.1332 

0.02623 

0.001882 

0.0000000 

6 

0.1336 

0.02631 

0.000332 

0.0000000 

7 

0.1343 

0.02652 

0.002509 

0.0000000 

8 

0.1340 

0.02608 

0.000279 

0.0000000 

9 

0.1341 

0.02659 

0.001415 

0.0000000 

10 

0.1326 

0.02628 

0.001764 

0.0000000 

11 

0.1343 

0.02641 

0.001312 

0.0002664 

12 

0.1339 

0.02650 

0.002210 

0.0000549 

13 

0.1344 

0.02617 

0.000710 

0.0000000 

14 

0.1332 

0.02537 

0.002150 

0.0000000 

15 

0.1334 

0.02670 

0.002611 

0.0000000 

16 

0.1328 

0.02683 

0.002629 

0.0000000 

17 

0.1334 

0.02662 

0.002360 

0.0001910 

18 

0.1336 

0.02686 

0.001415 

0.0000000 

19 

0.1341 

0.02619 

0.002322 

0.0000000 

20 

0.1341 

0.02636 

0.002085 

0.0000224 

21 

0.1337 

0.02680 

0.002059 

0.0000000 

22 

0.1342 

0.02615 

0.002300 

0.0002322 

23 

0.1337 

0.02464 

0.000000 

0.0000000 

24 

0.1343 

0.02566 

0.001393 

0.0000000 

25 

0.1342 

0.02677 

0.002103 

0.0002603 

max 

0.1344 

0.02686 

0.002629 

0.0002664 
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Figure  4.5:  The  angular  velocity  response  uj  =  Y,ty(u)c)  and  the  command  ujc  £ 
527r-5o(1350)  with  system  gain  G  =  10  are  shown  for  the  18  sequence  in  the  Monte 
Carlo  run. 
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Figure  4.6:  The  absolute  angular  velocity  error,  \tock  —D  2Tl^(cock)\,  is  shown  for  the 
18th  sequence  in  the  Monte  Carlo  run  of  qc  G  52^50(1350)  with  system  gain  G  —  10 
where  \uc  -  D-2lN(u>c)\  =  0.02686. 
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Figure  4.7:  The  quadrature  voltage  vq  is  shown  for  the  18th  sequence  in  the  Monte 
Carlo  run  of  u>c  €  52^.50(1350)  with  system  gain  G  =  10. 
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Figure  4.8:  The  absolute  quadrature  voltage  error,  \v'    —  vqk\,  is  shown  for  the  18 


th 


sequence  in  the  Monte  Carlo  run  of  ujc  £  527r.5o(1350)  w^h  system  gain  G  —  10  where 
\v'-va\  =  0.0993. 
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where  the  system  £  and  an  approximate  right  inverse  system  ^  are  illustrated  in 
Figure  4.3  and  Figure  4.4,  respectively. 

The  first  block  of  Figure  4.4  calculates  the  raw  quadrature  voltage  v'  for  a 
given  angular  velocity  command  uJck-  In  the  GDI  block,  the  voltage  quadrature 
command  vqmVk  is  determined  by  using  (4.22)  on  the  raw  quadrature  voltage  v'qk.  In 
Figure  4.3  (vcqk  —  vqinVk),  the  quadrature  voltage  vqk  is  the  output  of  the  DZ  block 
by  using  (4.21).  The  composition  of  (4.22)  and  (4.21)  yields 


vqk  =  DZ(CDl(v'qk)) 


It  follows  then  that  for  v'    =  ±^, 


{< 

if 

<  >  ch 

G-v'qk-l 

if 

h  <  <  <  ch 

0 

if 

G  —  uqk  —  G 

G.t^  +  1 

if 

G-l   <  <k  <  - 

1 
G 

I  <k 

if 

<<        G-l' 

(4.24; 


sup 
v'es(X) 


<?k 


Vq  ~  V1 


sup 
v'esm 


v'q-DZ(CDl(v'qk)) 


1 
G 


(4.25) 


For  G  =   10,   ^   =  0.1,  and  this  is  consistent  with  Figure  4.8  where  \v'  —  vq\   = 
0.09993  <  £. 

The  maximum  quadrature  voltage  error  occurs  when  the  magnitude  of  the  raw 
quadrature  voltage  \v'qk\  =  ^.  Since  there  is  a  feedback  loop  inside  an  approximate 
inverse  system  $,  the  quadrature  voltage  errors  are  not  accumulative  with  respect  to 
the  angular  velocity  errors.  Evaluating  the  right  hand  side  of  (4.23)  with  \v'\  =  ^ 
will  give  us  the  estimate  of  the  performance  bound  which  requires  two  integrations 
of  the  errors  as  follows: 


\Av 


'ik 


or 


|A»W+1| 


1     1 1 '  I    l 
J-a* 

GL 

60       Km 

2irGLJ 


AV 


60      Kr 


-At2. 


(4.26) 
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Using  the  optimization  set  uoc  €:  5,2tt-5o(1350),  and  for  G  =  10,  this  yields  "P(H)  = 
0.02688.  The  estimate  is  0.00002  higher  than  the  value  from  Table  4.2  of  0.02686. 
The  estimate  of  the  performance  bound  V(Yj)  is  directly  proportional  to  the  system 
gain  G.  Figure  4.6  is  a  scaling  of  Figure  4.8  by  a  factor  of  |^^jA£2.  In  this  example, 
we  have  found  a  closed  form  solution  for  an  estimate  of  the  performance  bound 
given  by  V(T,)  which  was  verified  by  simulation.  In  the  next  2  examples,  nonlinear 
programming  algorithms  will  be  used  to  calculate  "P(S). 

4.4     Aerodynamic  Model  Example 

This  section  will  provide  an  example  of  the  performance  bound  calculation  for 
the  longitudinal  control  of  an  aircraft.  Control  of  pitch  attitude  of  an  aircraft  can 
be  achieved  by  deflecting  all  or  portion  of  either  a  forward  or  aft  tail  surface.  The 
derivation  of  the  equations  of  motion  of  the  aircraft  with  the  assumption  that  the 
roll  rate,  yaw  rate  and  y-velocity  component  in  the  body-axes  frame  are  all  zero  is 
found  in  [4].  These  equations  are: 

u    =     —\X(u,w,q,6)  +  T-mgsm{0)]-wq  (4.27) 

m 

v     =     —[Z{u,w,q,6)  +  mgcos{0)]  +  uq  (4.28) 

m 

q    =    M(u,w,q,0,6) /lyy  (4.29) 

0    =    q  (4.30) 

where  u  and  w  are  the  body-axis  ar-velocity  component  and  the  body-axis  z- velocity 
component,  respectively.  The  body-axis  pitch  rate  and  the  Euler  pitch  angle  are 
given  by  q  and  9,  respectively.  The  longitudinal  control  input  or  elevator  deflection 
is  given  by  8.  X  and  Z  are  the  aerodynamic  forces  and  M  is  aerodynamic  moment. 
The  following  are  constants:  T,  thrust;  m,  mass;  g,  gravity;  and  \yy,  moment  of 
inertia.   The  aerodynamic  forces  and  moments  are  represented  by  the  means  of  the 
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aerodynamic  stability  coefficients  as  seen  below: 


q     =     arctan(u;/w) 


Q  = 

Q  = 

L  = 

D  = 

X  = 

z  = 

M  = 


(cos(a))   (wu  —  wit] 


u- 


-p  (u2  +  w2) 

{CL  +  CLaa  +  CLqq  +  CL66)QS 

{CD  +  CDaa)QS 

Lsin(a)  —  Dcos(a) 

—  L  cos(a)  —  D  sin(a) 

(Cm  +  Cmaa  +  Cmaa  +  Cmqq  +  Cms6)  QSc. 


(4.31) 
(4.32) 
(4.33) 


Here  a  and  a  are  the  angle  of  attack  and  the  angle  of  attack  rate,  respectively. 
Q  is  the  dynamic  pressure.  The  aerodynamic  moment,  M,  and  the  aerodynamic 
forces,  X  and  Z,  are  directly  proportional  to  t  he  dynamic  pressure.  L  and  D  are  the 
aerodynamic  forces  in  the  wind  axes  which  are  known  as  lift  and  drag  .  The  following 
are  constants:  5,  the  wing  reference  area  and  c,  the  wing  mean  aerodynamic  cord. 


Table  4.3:  The  Aerodynamic  Model  Parameters 

Parameter 

Value 

Parameter 

Value 

m 

1247.390  kg 

cD 

0.05 

T 

1600  N 

cDa 

0.33 

g 

9.81  m-s~2 

CL 

0.41 

lyy 

4066  kg-m2 

cLa 

4.44 

P 

1.225  kg-m~3 

cLq 

3.8-c/(2w0) 

s 

17.094  m2 

Cls 

0.355 

c 

1.737  m 

^ma 

-0.683 

u0 

53.8931  m-s"1 

Cm6 

-4.36-c/(2w0) 

v0 

-0.0931634  m-s"1 

Cmq 

-9.96-c/(2m0) 

Qo 

0  rad-s-1 

^m6 

-0.923 

0o 

0  rad 

6_i 

0.00128510  rad 

The  reference  geometry,  mass  characteristics,  and  aerodynamic  characteris- 
tics of  the  general  aviation  airplane:   NAVION  were  taken  from  Appendix  B  of  the 
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reference  [33].  The  values  of  the  parameters  are  shown  in  Table  4.3.  The  equa- 
tions of  motion  are  discretized  using  second  order  Runge-Kutta  integration.  Let 
Xk  =  [uk  Vk  qk  &k]'.  Equations  (4.27)  through  (4.30)  take  the  following  form: 


Zfc+i     =    f{xk)+g{xk)'6k 


qk   = 


0    0    10 


Xk 


(4.34) 
(4.35) 


where  Xk  €  3?4  and  /,  g  are  vector  fields  on  3?4. 

The  actuator  has  a  first-order  response  modelled  by  G(s)  =  10/(5  +  10)  which 
has  rate  constraints  of  15  deg/sec  and  position  constraints  of  15  deg.  The  system 
E:D  — *  5(3?)  :  6c  >-*  q,  where  V  C  5(3?),  is  the  composition  of  the  equations  of 
motion  with  the  actuator  dynamics  for  which  6c  is  the  elevator  deflection  command. 
Figure  4.9  shows  the  block  diagram  of  the  system  E. 

4.4.1     Simulation 

The  signal  of  interest  is  the  pitch  rate  signal  q.  It  is  assumed  that  all  signals 
of  interest  start  at  rest  (i.e.,  qn  =  0  for  all  n  <  no).  The  absolute  bound  of  \q\  <  1.15 
deg/sec  was  imposed  so  that  after  10  seconds  of  flight,  the  aircraft  will  always  be 
within  its  parameters  of  good  flight  conditions.  Typical  pitch  rate  commands  are 
generated  by  the  longitudinal  autopilot.  The  altitude  hold  mode  and  the  mach  hold 
mode  have  time  constants  of  approximately  100  milli-seconds  ([4]).  This  corresponds 
to  a  restriction  on  bandwidth  of  q  to  be  approximately  10  rad/sec. 
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Rate  Position 

Constraints         Constraints 
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z~'         < 


xk+l  =  f(x  k)  +  g(x  k)  5k 
q'k  =  [0  0  1  0]  x   . 


jflfc+1 


Figure  4.9:  The  block  diagram  of  the  system  E:X>  — >  5(3?)  :  6c  *—>■  q  for  the  Aerody- 
namic Model  Example. 
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The  block  diagram  of  an  approximate  right  inverse  system  $:  SnLP(a)  — ♦ 
5(3?)  :  qc  |-+  <^inv  is  shown  in  Figure  4.10.  The  first  block  calculates  the  raw  elevator 
deflection  8'k  for  a  given  pitch  rate  command  qck.  In  the  limiting  logic  block,  the 
elevator  deflection  command  <!J;nv  is  determined  by  first  using  the  constraints  on  the 
raw  elevator  deflection  8k  and  then  using  the  inverse  actuator  dynamics.  The  elevator 
deflection  command  8mv  is  then  used  to  update  the  states  xinv  of  the  system  $. 

One  step  delay  of  8c  is  required  to  reach  q.  Therefore,  £(£)  >  1.  From  Fig- 
ure 4.10  it  follows  that  £($)  >  0.  Theorem  3.1.2  applied  to  the  composite  system 
E#:SfiLP(cr)  -J-  5(3?)  :  qc  >->  q  yields  £(E#)  >  1.  For  the  considered  flight  condi- 
tions, g(xk)  ^  0  in  (4.34)  for  all  8c  €  V  and  a  given  x0.  The  design  integer  is  taken 
as  A  =  1  for  an  approximate  right  inverse  system  $. 

For  a  simulation  of  10  seconds,  a  sampling  period  of  40  milli-seconds  yields 
a  sequence  of  250  points.  The  pre-filtered  binary  random  sequences  have  only  the 
values  ±1.15  deg/sec  because  we  are  only  interested  in  the  extreme  points  of  5(1.15). 
Table  4.4  shows  the  Monte  Carlo  results  for  3  different  cutoff  frequencies  Q,lp'-  10, 
15,  and  20  rad/sec,  where  the  error  \qc  —  D~l T,ty(qc)\  is  evaluated  for  each  of 
the  sequences.  For  practical  purposes,  as  can  be  seen  from  Table  4.4,  the  system 
#:  5io(1.15)  — >  5(3?)  is  a  right  inverse  of  the  system  S. 
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c  I 
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Figure  4.10:  The  block  diagram  of  an  approximate  right  inverse  system  $ :  Sqlp  (a) 
5(3?)  :  qc  >-*  8\nv  for  the  Aerodynamic  Model  Example. 
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Table  4.4:  Family  of  Sequences,  SnLF(1.15),  Comparison 


Sequence 

Approximate  Ri 

ght  Inverse 

Error  \qc  -  D-xJM(qc)\ 

Qlp  —  10  rad/sec 

VtLP  =  15 

rad/sec 

QLp  =  20  rad/sec 

1 

4.996E-16 

0.9075 

1.353 

2 

4.163E-16 

1.7650 

2.265 

3 

3.331E-16 

0.4226 

1.553 

4 

2.776E-16 

0.3620 

1.121 

5 

3.608E-16 

0.5654 

1.592 

6 

3.331E-16 

0.4446 

1.356 

7 

3.331E-16 

1.4340 

2.153 

8 

2.498E-16 

0.9444 

1.788 

9 

3.331E-16 

0.6539 

1.896 

10 

4.441E-16 

1.2530 

1.931 

11 

2.220E-16 

1.2270 

2.121 

12 

5.551E-16 

0.9812 

1.371 

13 

4.996E-16 

0.7821 

1.471 

14 

5.551E-16 

0.7177 

1.488 

15 

4.441E-16 

0.6425 

1.351 

16 

4.441E-16 

0.7448 

1.602 

17 

5.551E-16 

0.4841 

1.454 

18 

4.996E-16 

0.4835 

1.586 

19 

4.441E-16 

0.6658 

1.585 

20 

3.331E-16 

0.7861 

1.675 

21 

3.608E-16 

0.5702 

1.544 

22 

3.331E-16 

0.6071 

1.809 

23 

3.886E-16 

0.6047 

1.635 

24 

3.053E-16 

0.4456 

1.797 

25 

4.441E-16 

0.9205 

1.413 

max 

5.551E-16 

1.7650 

2.265 
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Figure  4.11:    The  shifted  pitch  rate  response  q  =  D  1T,1$(qc)  and  the  command 
qc  £  5,i5(1.15)  are  shown  for  the  2nd  sequence  in  the  Monte  Carlo  run. 
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Figure  4.12:  The  absolute  pitch  rate  error,  \qck  —  D  1E$(qck)\,  is  shown  for  the  2nd 
sequence  in  the  Monte  Carlo  run  of  qc  £  5is(1.15)  where  \qc  —  D~1T,^(qc)\  =  1.765. 
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Figure  4.13:    The  elevator  deflection  6  is  shown  for  the  2nd  sequence  in  the  Monte 
Carlo  run  of  qc  G  5i5(1.15). 
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Figure  4.14:  The  elevator  deflection  rate  8  is  shown  for  the  the  2nd  sequence  in  the 
Monte  Carlo  run  of  qc  €  5is(1.15). 
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From  Table  4.4,  the  maximum  value  of  \qc  —  D~1H1$!(qc)\  for  f^p  =  15 
rad/sec  is  1.765.  This  corresponds  to  the  2nd  sequence  in  the  Monte  Carlo  run.  In 
the  following  subsection  we  will  calculate  the  estimate  of  the  performance  bound  "P(S) 
over  the  optimization  set  615(1. 15).  A  first  estimate  is  V{Ti)  >  1.765.  In  Figure  4.11, 
the  2nd  sequence  of  qc,  the  filtered  square  wave,  is  plotted  with  the  corresponding 
shifted  sequence  T,^(qc)  which  is  shifted  one  step  to  the  left  by  40  milli-seconds.  The 
absolute  pitch  rate  error  sequence  is  plotted  in  Figure  4.12.  The  elevator  deflection  6 
and  elevator  deflection  rate  6  are  plotted  in  Figure  4.13  and  Figure  4.14,  respectively. 
As  can  be  seen  from  Figure  4.12,  the  maximum  error  occurs  at  8.0  seconds.  Examining 
the  plots  between  7  and  8  seconds  reveals  that  the  pitch  rate  response  is  not  keeping 
up  with  the  pitch  rate  command.  The  elevator  deflection  is  never  position  constrained 
but  is  rate  constrained  from  7.44  to  8.12  seconds  for  a  total  of  680  milli-seconds.  As 
seen  from  this  example,  when  we  increase  the  bandwidth  of  the  pitch  rate  command 
from  10  to  15  rad/sec,  the  rate  constraints  become  effective  and  thus  the  measure  of 
right  singularity  is  increased. 

4.4.2     The  Estimate  of  the  Performance  Bound  V(T,) 

In  this  section,  we  will  find  an  estimate  of  the  performance  bound  "P(S)  for 

the  longitudinal  control  problem  of  an  aircraft  by  using  the  optimization  definition 

as  follows: 

V(Y)=        sup        qc-D'^^iqc)  (4.36) 

<?ces15(i.i5) 

where  the  system  E  and  an  approximate  right  inverse  system  ^  are  illustrated  in 
Figure  4.9  and  Figure  4.10,  respectively.  The  discussion  following  Theorem  4.2.1 
in  Section  4.2  suggests  that  the  maximum  of  the  right  hand  side  of  (4.36)  is  found 
by  evaluating  all  extreme  points  and  then  finding  the  maximum.  In  Section  4.4.1, 
25  random  extreme  points  out  of  possible  2250  (1.8903  x  1075)  extreme  points  were 
evaluated  with  the  conclusion  P(E)  >  1.765. 
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q    =1.15 


<?    =-1.15 


hNp, 


-Np. 


Figure  4.15:  Signal  Construction's  Nomenclature. 

It  is  time  to  take  advantage  of  the  traits  of  the  system.  The  first  block  of 
Figure  4.10  calculates  the  raw  elevator  deflection  S'k  for  a  given  pitch  rate  command 
qck.  In  mathematical  terms  it  is  solving  for  S'k  in  the  following  equations: 

[<lck  ~qk-\) 


Ml     =    I 


yy 


At 


M[      =       {Cm  +  CmaQk-l   +  Cmqqk-i  + 

Cmaa(uk-\,Wk-\,8'k)  +  Cm^'k)Qk-\Sc. 


(4.37) 


(4.38) 


As  can  be  seen  from  (4.38),  the  raw  aerodynamic  moment  M'k  is  affme  with  respect 
to  the  dynamic  pressure  Qk-i-  It  turns  out  that  it  is  also  affme  with  respect  to  the 
raw  elevator  deflection  S'k.  Thus,  S'k  is  inversely  proportional  to  Qk-\  and  decreasing 
dynamic  pressure  causes  the  elevator  deflection  to  increase.  This  increases  the  likeli- 
hood for  the  elevator  deflection  to  be  rate  constrained.  Therefore,  the  optimization 
region  of  (4.36)  can  be  reduced  to  a  subset  of  5i5(1.15)  which  contains  a  lower  than 
initial  dynamic  pressure  after  a  period  of  flight.  A  pitch-up  command,  say  a  com- 
mand qc  =  1.15  deg/sec  for  at  least  2  seconds  after  initialization,  will  decrease  the 
dynamic  pressure. 


53 

Two  problems  of  maximization  of  convex  functions  are: 

1.  Once  a  local  maximum  has  been  found.    There  is  no  local  information  to  tell 
you  how  to  proceed  to  a  higher  local  maximum. 

2.  There  is  no  local  criterion  for  deciding  whether  a  given  local  maximum  is  really 
the  global  maximum. 

These  features  are  also  applicable  to  the  maximization  of  non-convex  function  since  its 
maximum  for  simulation  purposes  will  be  an  extreme  point.  The  signal  construction 
and  the  grid  search  method  will  be  techniques  used  to  find  a  higher  local  maximum 
to  get  around  the  first  problems.  These  two  methods  are  explained  in  the  following 
paragraphs.  After  an  extensive  search  for  the  signal  that  maximizes  (4.36)  and  an 
appropriate  exiting  condition,  we  can  say  that  we  have  qualitatively  found  the  global 
maximum.  This  alleviates  the  second  problem. 

The  signal  construction  will  always  be  done  in  the  pre-filtered  space  5(1.15) 
with  out  any  loss  of  generality  since  the  pre-filtered  space  5(1.15)  is  isomorphic  to 
the  filtered  space  5is(1.15).  The  baseline  signal  will  be  qc  =  1.15  deg/sec  for  the 
entire  simulation  of  10  seconds.  The  construction  of  the  signals  begins  with  a  pulse 
qc  —  —1.15  for  Np\  steps  and  starting  at  N\  steps.  An  optimization  takes  place 
over  Ni  and  Npx.  The  signal  construction's  nomenclature  is  shown  in  Figure  4.15.  A 
second  pulse  is  constructed  with  Nu2  steps  at  qc  —  1.15  and  Np2  steps  at  qc  =  —1.15. 
The  optimization  takes  place  over  N2,  Nu2,  and  Np2  with  Npi  fixed.  This  procedures 
continues  until  the  (i  +  1  )-th  pulse  optimization  does  not  produce  higher  maxima. 

A  signal  with  i  pulses  has  2z  parameters,  (iV,-,  Npi, . . . ,  Npi,  Nu2, . . . ,  Nul).  In 
the  search  grid  method  the  starting  parameter  N{  (see  Figure  4.15)  of  the  pulse  train 
is  still  optimized  over  its  range.  The  other  2i  —  1  parameters  (Npi, . . . ,  Npi,  Nu2, . . . , 
Nui)  are  varied  over  a  selected  range.  The  selected  range  is  based  on  the  available 
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Figure  4.16:    The  shifted  pitch  rate  response  q  =  D  lY,ty(qc)  and  the  command 
qc  G  5i5(1.15)  are  shown  for  the  first  pulse  optimization  {N\  =  242,  Npi  =  4). 
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Figure  4.17:  The  absolute  pitch  rate  error,  \qck-  D  lT,ty{qck)\,  is  shown  for  the  first 
pulse  optimization  (iVi  =  242,  Npl  =  4)  of  qc  €  Sis(1.15)  where  \qc  -  D^IM  (qc)\  = 
1.646. 


55 


4  5  6 

Time  (seconds) 


10 


Figure  4.18:    The  shifted  pitch  rate  response  q  =  D   1Tl^(qc)  and  the  command 
qc  G  5i5(1.15)  are  shown  for  the  second  pulse  optimization  (N2  =  222,  Nu2  =  7, 

Np2  =  7). 
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Figure  4.19:  The  absolute  pitch  rate  error,  \qck  —  D  1E^(^cfe)|5  is  shown  for  the 
second  pulse  optimization  (JV2  =  222,  Nu2  =  7,  Np2  =  7)  of  qc  G  5*1 5 ( 1 . 1 5 )  where 
\qc-D-lm(qc)\  =  2.713. 
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Figure  4.20:  The  shifted  pitch  rate  response  q  =  D  1Hty(qc)  and  the  command 
qc  £  -Si 5 ( 1 . 1 5 )  are  shown  for  the  third  pulse  optimization  (N3  =  195,  Nu3  =  8, 
NP3  =  7). 
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Figure  4.21:  The  absolute  pitch  rate  error,  \qck  —  D  1T,ty(qck)\,  is  shown  for  the 
third  pulse  optimization  {N3  =  195,  Nu3  =  8,  Np3  -  7)  of  qc  6  5i5(1.15)  where 
\qe-D-lVW(qc)\  =  2.856. 
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Figure  4.22:  The  shifted  pitch  rate  response  q  —  D  lTi^!{qc)  and  the  command  qc  G 
5i5(1.15)  are  shown  for  the  search  grid  method  optimization  (iV3  =  215,  Npi  =  2, 
Nu2  =  7,  Np2  =  5,  Nu3  =  7,  Np3  =  8). 
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Figure  4.23:  The  absolute  pitch  rate  error,  \qck 
search  grid  method  optimization  (N3  =  215,  ./Vpi 
Np3  =  8)  of  qc  G  Si5(1.15)  where  P(E)  =  2.908. 


D   1T,]$(qck)\,  is  shown  for  the 
2,  Nu2  =  7,  Np2  =  5,  Nu3  =  7, 
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Figure  4.24:  The  elevator  deflection  8  is  shown  for  the  search  grid  method  optimiza- 
tion {N3  =  215,  NPl  =  2,  Nu2  =  7,  Np2  =  5,  Nu3  =  7,  Np3  =  8)  of  qc  €  515(1.15). 
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Figure  4.25:  The  elevator  deflection  rate  8  is  shown  for  the  the  search  grid  method 
optimization  (JV3  =  215,  Npl  =  2,  Nu2  =  7,  Np2  =  5,  Nu3  =  7,  Np3  =  8)  of 
qc  e  5W(1.15). 
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computing  power.  All  points  in  the  selected  range  are  evaluated  and  then  the  max- 
imum is  found.  The  maximum  of  the  search  grid  method  is  qualitatively  the  global 
maximum  provided  that  the  selected  range  is  suitable. 

The  results  of  the  nonlinear  programming  algorithm,  which  is  the  integration 
of  the  signal  construction  and  the  grid  search  method,  are  shown  in  Table  4.5  and 
Figures  4.16-4.25.  Figure  4.16  and  Figure  4.17  represent  the  first  pulse  optimization 
with  plots  of  the  sequences  qc,  D~1Y,ty(qc)  and  \qck  —  D~1Y,ty(qck)\-  Figure  4.18 
and  Figure  4.19  represent  the  second  pulse  optimization  with  plots  of  the  sequences 
qc,  D_1E*(9c)  and  \qCk  -  D~lYNl(qch)\.  Figure  4.20  and  Figure  4.21  represent 
the  third  pulse  optimization  with  plots  of  the  sequences  qc,  D~1Y,ty(qc)  and  \qck  — 
D~^{qck)\. 

The  third  pulse  optimization  parameters  are  {Np\,  Nu2,  Np2,  Nv.3,  Np3)  = 
(4,7,7,8,7).  The  search  grid  was  taken  to  be  (2:6,5:9,5:9,6:10,5:9)  which  yields 
3125  grid  points  that  have  to  be  optimized  over  the  variable  A3.  The  result  of  the 
search  grid  method  and  the  estimate  of  the  performance  bound  are: 

(AT3,  NPl,Nu2,  Np2,  Nu3,  Np3)    =    (215,  2,  7,  5,  7,  8) 

P(E)    =    2.908. 

The  search  grid  sequence  of  qc  €  5i5(1.15)  is  examined.  In  Figure  4.22,  the  se- 
quence qc,  the  filtered  square  wave,  is  plotted  with  the  corresponding  shifted  se- 
quence Tity(qc)  which  is  shifted  one  step  to  the  left  by  40  milli-seconds.  The  absolute 
pitch  rate  error  sequence  is  plotted  in  Figure  4.23.  The  elevator  deflection  8  and 
elevator  deflection  rate  6  are  plotted  in  Figure  4.24  and  Figure  4.25,  respectively.  As 
can  be  seen  from  Figure  4.23,  the  maximum  error  occurs  at  9.88  seconds.  Examin- 
ing the  plots  between  8.6  and  10  seconds  reveals  that  the  pitch  rate  response  does 
not  keep  up  with  the  pitch  rate  command.  The  elevator  deflection  is  never  position 
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constrained  but  is  rate  constrained  for  a  total  1.2  seconds  out  of  the  last  1.4  seconds 
of  simulation. 

Table  4.5:  The  Estimate  of  the  Performance  Bound  'P(E)  for  the  Aerodynamic  Model 
Example 


Optimization 

Method 

Parameters 

(TV,,  NPl ,  Nu2,  Np2,  Nu3,  Np3) 

\qc  -  D-lW{qc)\ 
for  qc  e  5i5(1.15) 

First  Pulse 
Second  Pulse 
Third  Pulse 
Search  Grid 

(242,4) 
(222,4,7,7) 
(195,4,7,7,8,7) 
(215,2,7,5,7,8) 

1.646 
2.713 
2.856 
2.908 

P(E)  =2.908 

4.5     Multivariate  Process  Control  Model  Example 

This  section  is  concerned  with  multivariable  control  of  a  plant  consisting  of  a 
cylindrical  tank  containing  a  water  column  pressurized  by  air.  This  example  comes 
from  [10,  11]  wherein  the  authors  developed  a  multivariable  predictive  control  strategy 
with  a  variable  receding  horizon. 

The  process  control  model  consists  of  a  cylindrical  tank  which  contains  a  water 
column  pressurized  by  air.  Figure  4.26  is  a  simplified  process  diagram.  A  centrifugal 
pump  feeds  water  into  the  tank  through  the  pneumatic  valve  1,  and  a  compressor 
feeds  air  through  the  pneumatic  valve  2.  Liquid  can  drain  from  the  tank  through 
the  manual  valve  3  while  air  can  exit  through  the  relief  valve  4.  Measurements  of 
the  process  output — the  height  of  the  liquid  column  y\  and  the  gauge  pressure  of 
the  enclosed  air  y2 — are  obtained  by  means  of  pressure  transducers.  Valves  1  and  2 
are  manipulated  through  electro-pneumatic  transducers  by  voltage  signals  U\  and  u2. 
The  height  of  the  tank  is  105  cm  and  the  diameter  is  7.5  cm. 

The  dynamics  of  the  process  control  model  can  be  described  in  terms  of  four 
states.  Let  xx  denote  the  liquid-level  height  [m],  x2  the  air  pressure  inside  the  column 
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Figure  4.26:  The  Multivariable  Process  Control  Model  with  inputs:  valve  signals  ut 
and  w2;  with  outputs:  liquid  level  y\  and  the  air  pressure  y2. 

[bar  gauge],  £3  the  opening  of  water  valve  1  (lift)  [dimensionless],  and  x4  the  opening 
of  air  valve  2  (lift)  [dimensionless].  The  equations  of  motion  of  the  process  control 
model  are  given  by: 

1 


X!  = 

X2  = 

?3  = 

X4  — 


(9//  -  qid) 


i  +  z2 


S(L-Xl) 

-  [gm  -  x3) 

—  (ggu2  -  x4) 


[{qif  -qid)  +  {q9j  -  qgdj\ 


(4.39) 
(4.40) 
(4.41) 
(4.42) 


where  L  is  the  height  of  the  column  and  S  is  the  cross-sectional  area  of  the  column. 
The  volumetric  flow  rates  of  the  liquid  feed  and  liquid  discharge  streams  are  repre- 
sented by  the  variables  qij  and  qid,  respectively;  and  the  gas  feed  and  gas  discharge 
rates  by  the  variables  qgj  and  qgd,  respectively. 

The  liquid  and  air  flow  rates  are  modeled  using  a  steady-state  relationships 
appropriate  for  the  square-root  and  equal-percentage  characteristics  of  the  valves 
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used  in  the  process  control  model: 


qif  =  x3KifyJpif  -  x2  -  piG{xx  +  AL') 

qid  =  xldKidsJx2  +  piG(x1  +  AL") 

qgf  =  x4Kgf{pgI  -  x2)0-65 

qgd  -  Kgd\/x~2 


(4.43) 
(4.44) 
(4.45) 
(4.46) 


where  the  constants  K{3  represent  the  valve  coefficients,  G  is  the  acceleration  of 
gravity,  pi  is  density  of  the  liquid,  pgf  is  the  inlet  pressure  of  the  air,  pij  is  the  inlet 
pressure  of  the  liquid,  AL'  is  the  height  difference  between  valve  1  and  the  bottom  of 
the  column,  and  AL"  is  the  height  difference  between  valve  3  and  the  bottom  of  the 
column.  Specific  values  and  definitions  of  all  the  multivariate  process  control  model 
parameters  are  given  in  Table  4.6.  The  states  and  inputs  must  satisfy  the  constraints 


"  0  1 

"  XX   ' 

"   Xi    " 

0 
0 

< 

X2 
X3 

< 

x2 

1 

1 

0 
0 

< 

u2 

< 

u2 

0 

£4 

1 

(4.47) 


where  u1  =  u2  —  10[V]  (the  maximum  signal  voltage)  and  X\  and  x2  are  bounds 
established  by  the  user.  The  bounds  of  the  output  vector  are  identical  to  those  of 
the  states  Xi  and  x2  because  y\  =  x^  and  y2  =  x2. 


Table  4.6: 

The 

'  Multivariate 

Process  Control  Model  Parameters 

Parameter 

Value 

Parameter 

Value 

S 

0.0043  m2 

Ku 

4.0  x  10-4  m3-s-1-v/bar 

L 

1.05  m 

Kit 

6.9  x  10~4  m^s-^-v/bar 

AL' 

0.16  m 

Kgd 

9.8  x  10-4  m3-s-1-v/bar 

AL" 

0.40  m 

K9f 

3.1  x  10~4  m3-s-1Vbar 

Pif 

0.5  bar 

n 

1.5  s 

Pi 

1.0 

x  103  kg-m-3 

9i 

0.1  v-1 

Pgf 

0.4  bar 

T9 

9g 

0.3  s 

0.1  v-1 
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The  equations  of  motion  are  discretized  using  first  order  Euler  integration. 
Let  Xk  =  [xik  x2k  x3k  x4J'.  Equations  (4.39)  through  (4.42)  take  the  following  form: 


Xk+i    =   f{xk)+    g\{xk)   gi{xk) 


u2k 


10    0    0 
0    10    0 


Xk 


(4.48) 
(4.49) 


where  Xk  €  9ft4,  and  /,  gi,  g2  are  vector  fields  on  3ft4.  The  system  T,:V  — ►  5(3ft2)  : 
uc  i— ►  ?/  of  interest,  where  V  C  5(3ft2),  is  the  composition  of  the  equations  of  motion 
with  the  input/states  constraints  for  which  uq  is  the  voltage  command  signals  U\  and 
u2.  Figure  4.27  shows  the  block  diagram  of  the  system  E. 
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Figure  4.27:    The  block  diagram  of  the  system  E:P  — >  5(3ft2)  :    uc  *— >■  y  for  the 
Multivariable  Process  Control  Model  Example. 


4.5.1     Simulation 

The  signals  of  interest  are  the  liquid-level  height  y\  and  the  gauge  pressure 
t/2  of  the  enclosed  air.  The  bounds  of  (4.47)  have  been  set  to  X\  =  0.75  m  and 
x2  =  0.30  bar.  It  is  assumed  that  all  signals  of  interest  start  at  their  midpoints  (i.e., 
j/in  =  0.375  m  and  y2n  =  0.15  bar  for  all  n  <  n0).  From  [11],  the  time  constant  of 
the  liquid-level  y\  (under  zero  gauge-pressure  of  air)  varies  from  50-70  s  and  the  time 
constant  of  the  air  pressure  y2  is  1.3  s.  The  bandwidth  of  y\  was  set  at  0.3  rad/sec 
to  stress  an  approximate  right  inverse  system  $  while  the  bandwidth  of  y2  was  set 
at  1/1.3  rad/sec. 
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The  set  T>s  C  £(3ft2)  is  defined  with  initial  conditions  j/i,  =  0.375  and  y2|  = 
0.15  for  all  i  <  0,  and  with  \yu  -  0.375|  <  0.357  and  \y2i  -  0.15|  <  0.1428,  for  all 
i  >  0.  The  domain  £\p  C  5(3ft2)  of  an  approximate  right  inverse  system  $  takes  the 
sequences  of  "D5's  components  1  and  2  that  are  passed  through  a  second  order  discrete- 
time  low  pass  filter  with  cutoff  frequencies  0.3  and  1/1.3  rad/sec,  respectively.  The 
outputs  of  the  filters  1  and  2  are  then  clipped  by  [0.0, 0.75]  and  [0.0, 0.30],  respectively. 
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Figure  4.28:  The  block  diagram  of  an  approximate  right  inverse  system  ty:V$  — * 
S($.2)  :  yc  •— >  w;nv  with  design  integer  A  =  2  for  the  Multivariate  Process  Control 
Model  Example. 
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Figure  4.29:  The  block  diagram  of  an  approximate  right  inverse  system  ty:T>y  — > 
5(3?2)  :  yc  !— >  "inv  with  design  integer  A  =  3  for  the  Multivariate  Process  Control 
Model  Example. 


The  block  diagrams  of  approximate  right  inverse  systems  ty:Vy   — >  5(3ft2) 
:    Vc  |— *  «inv  with  design  integers  A  =  2  and  A  =  3  are  shown  in  Figure  4.28  and 
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Figure  4.29,  respectively.  Two  integrations  of  u\  are  required  to  reach  y\  and  sim- 
ilarly for  u2  and  y2.  Two  integrations  of  U\  are  also  required  to  reach  y2  but  three 
integrations  of  u2  are  required  to  reach  y\. 

For  the  design  integer  A  =  2  case,  the  system  ty  is  causal.  The  u\  optimization 
block  in  Figure  4.28  optimizes  U\k  to  give  a  minimum  value  of  \yck  —  D~2T,ty(yck)\  each 
step.  For  the  design  integer  A  =  3  case,  the  system  $  is  strictly  causal.  The  U\  and  u2 
optimization  block  in  Figure  4.29  optimizes  U\k  and  u2k  to  give  a  minimum  value  of 
\yck  —D~3 Y,ty(yck) |  each  step.  The  procedure  involves  undertaking  a  U2k  optimization 
if  the  Uik  optimization  fails  so  as  to  reduce  the  error  \yck  —  D~3T,ty(yck)\-  The  signal 
u2k  is  optimized  to  give  a  minimum  value  of  \yck+l  —  D ~3T,ty(yck+1  )|  provided  the 
result  does  not  augment  the  error  \yck  —  D~3Y,ty(yck)\. 

For  the  simulation  of  60  seconds,  a  sampling  period  of  150  milli-seconds  yields 
a  two-dimensional  sequence  of  400  points.  The  pre-filtered  binary  random  sequences 
{ywk}  only  have  values  of  0.018,  0.0732  m  for  A;  >  1,  and  the  pre-filtered  binary 
random  sequences  {y2Ck}  only  have  values  of  0.0072,  0.2928  bar  for  k  >  1  because  we 
are  only  interested  in  extreme  points  of  Vs.  Table  4.7  shows  the  Monte  Carlo  results 
for  2  different  design  integers  A  of  the  system  $:  A  =  2  and  A  =  3  where  the  error 
\yc  —  D~x^,^(yc)\  is  evaluated  for  each  of  the  sequences.  Also  included  in  Table  4.7 
is  the  average  number  of  floating  point  operations  (flops)  per  Monte  Carlo  sequence 
run  for  the  Ui  optimization  block  and  for  the  U\  and  u2  optimization  block.  On 
the  average  the  Ui  and  u2  optimization  subroutine  requires  3x  greater  the  number 
of  flops  as  the  Ui  optimization  subroutine,  but  a  single  execution  of  the  u^  and  u2 
subroutine  can  require  4+  x  greater  the  number  of  flops  than  the  u\  subroutine. 

4.5.2     The  Estimate  of  the  Performance  Bound  ViT.) 

In  this  section,  we  will  find  an  estimate  of  the  performance  bound  V(T>)  for 
the  multivariate  process  control  problem  of  a  regulator  that  regulates  the  liquid  level 
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Table  4.7:  Family  of  Sequences,  P<i»,  Comparison 

Sequence 

Error  \yc   -  D-Xm(yc)\ 

Average  Number  of  Flops 

A  =  2 

A  =  3 

A  =  2     A  =  3 

1 

0.02481 

0.02190 

137.0      367.8 

2 

0.03742 

0.03138 

140.3      449.3 

3 

0.01584 

0.01506 

142.5      396.0 

4 

0.05440 

0.04699 

150.0      432.7 

5 

0.02408 

0.02255 

135.5      447.5 

6 

0.02064 

0.01732 

137.3      441.0 

7 

0.02147 

0.01917 

136.1      379.4 

8 

0.02613 

0.02389 

140.9      426.5 

9 

0.01185 

0.01124 

135.5      401.5 

10 

0.00601 

0.00589 

139.7      374.2 

11 

0.02080 

0.02089 

144.9      459.0 

12 

0.03080 

0.02840 

137.0      439.1 

13 

0.01270 

0.01113 

139.1      404.9 

14 

0.02746 

0.02446 

135.5      456.5 

15 

0.03962 

0.03634 

144.0      513.7 

16 

0.02189 

0.01942 

140.3      413.9 

17 

0.01852 

0.01713 

134.6      454.2 

18 

0.01745 

0.01819 

136.1      467.9 

19 

0.01757 

0.01527 

138.5      441.9 

20 

0.01637 

0.01393 

143.1      399.8 

21 

0.03920 

0.03406 

144.6      435.2 

22 

0.03529 

0.02838 

138.2      418.2 

23 

0.01935 

0.01764 

140.9      464.1 

24 

0.01943 

0.01733 

138.5      420.2 

25 

0.01594 

0.01466 

139.4      393.5 

max 

0.0544 

0.04699 

150.0      513.7 

min 

134.6      367.8 

mean 

139.6      427.9 

st  dev 

3.703      33.90 
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Figure  4.30:   For  A  =  2,  the  liquid  height  response  y\  =  D  2E^(t/ic)  and  the  com- 
mand y\c  £  T>q  are  shown  for  the  square  wave  optimization  (A^  =  104,  N\2  =  103). 
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Figure  4.31:   For  A  =  3,  the  liquid  height  response  yx  =  D  3E$(?/ic)  and  the  com- 
mand y\c  £  £>#  are  shown  for  the  square  wave  optimization  {N^  =  78,  iVl2  =  73). 
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Figure  4.32:  For  A  =  2,  the  air  pressure  response  y2  =  D  2E^(j/2c')  and  the  command 
V2C  €  T>y  are  shown  for  the  square  wave  optimization  {N\1  =  104,  N\2  =  103). 


0.3 


—  0.25 


™ 

E 
E 
o 
O 

T3 

C 
CO 


0.2- 


0.15 


£    0.1 


CD 


0.05 


r  "■^  1 

1                               1                               1 

1 

/    1 

—  y2 

'                 /  \ 

'/           ' 
'/ 

I                       1 
J                        l' 
J                        1 

\ 

-  -   y2C 

- 

1 

I 

1 

1 

1 

1 

1 

\ 
1 
1 
\ 

\ 

,   1                          1                          1 

1 

10 


20 


30 
Time  (seconds) 


40 


50 


60 


Figure  4.33:  For  A  =  3,  the  air  pressure  response  y2  =  D  3E^(y2c)  and  the  command 
J/2C  £  2\  are  shown  for  the  square  wave  optimization  (TVij  =  78,  N\2  =  73). 
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Figure  4.34:  For  A  =  2,  the  absolute  errors,  \\)\ck  —  D  2H^(yick)\  and  \y2Ck  — 
D~2Tity(y2ck)\,  are  shown  for  the  square  wave  optimization  (Ni1  —  104,  7Vi2  =  103) 
of  yc  €  X>g>  where  V{T.)  =  0.5983  and  \y2C  -  D-2Y1^{y2C)\  =  0.0059. 
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Figure  4.35:  For  A  =  3,  the  absolute  errors,  \y\ck  ~  D  3T,ty(yick)\  and  \y2ck  ~ 
D~3T,$(y2ck)\1  are  shown  for  the  square  wave  optimization  {Ni1  =  78,  7Vi2  =  73)  of 
yc  €  £>*  where  P(S)  =  0.4544  and  \y2C  -  £>-3E#(y2C7)|  =  0.2495. 
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Figure  4.36:  For  A  =  2,  the  valve  commands  «i  and  w2  are  shown  for  the  square  wave 
optimization  (Nu  =  104,  Nl2  =  103)  of  yc  €  £>*. 
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Figure  4.37:  For  A  =  3,  the  valve  commands  U\  and  u2  are  shown  for  the  square  wave 
optimization  (TVi,  =  78,  N\2  =  73)  of  yc  £  T>y. 
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in  a  pressurized  tank  by  using  the  optimization  definition  as  follows: 


P(E)  =    sup 

yc£Vy 


yc-D-xW(yc 


(4.50) 


where  the  system  £  and  an  approximate  right  inverse  system  $  with  design  integers 
A:  2  and  3  are  illustrated  in  Figure  4.27,  Figure  4.28,  and  Figure  4.29,  respectively. 
In  Section  4.4.2,  a  nonlinear  programming  algorithm  was  developed  to  calcu- 
late the  estimate  of  the  performance  bound  "P(S).  It  exploits  the  traits  of  the  system. 
Due  to  the  sluggishness  of  the  liquid  height  response,  only  one  transition  was  needed 
from  high  to  low  in  both  components  of  the  signal  yc-  The  transition  is  at  iV^  steps 
for  yic  and  7Vi2  steps  for  y2c-  The  results  of  the  nonlinear  programming  algorithm 
are  shown  in  Table  4.8  and  Figures  4.30-4.37. 


Table  4.8:  The  Estimate  of  the  Performance  Bound  "P(E)  for  the  Multi variable  Pro- 
cess Control  Model  Example 


Design  Integer 

Parameters 
(Nu,Nl2) 

\yc 

-D-^(yc)\ 
for  yc  £  Vy 

A  =  2                   (104,103) 
A  =  3                    (78,73) 

P(E)  =  0.5983 
P(E)  =  0.4544 

Figure  4.30  and  Figure  4.31  represent  the  optimal  liquid  height  response  and 
command  for  A  =  2  and  A  =  3,  respectively.  The  U\  and  u2  subroutine  improves 
the  liquid  height  response  for  A  =  3,  at  the  expense  of  the  air  pressure  response. 
Figure  4.32  and  Figure  4.33  represent  the  optimal  air  pressure  response  and  command 
for  A  =  2  and  A  =  3,  respectively.  The  u\  and  u2  subroutine  increases  the  air  pressure 
to  help  drain  the  tank  while  the  iti  subroutine  is  operating  with  near  zero  gauge 
pressure  on  y\  which  has  a  time  constant  between  50-70  s. 

The  optimal  absolute  error  sequences  are  plotted  in  Figure  4.34  and  Fig- 
ure 4.35  for  A  =  2  and  A  =  3,  respectively.  In  the  u\  h  u2  subroutine,  the  signal  u2k 
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is  optimized  to  give  a  minimum  value  of  \yck+1  —  D~3Y,ty(yck+l  )|  provided  the  error 
\yck  —  D~3Ti^{yck)\  does  not  augment.  In  the  process,  the  error  \y2ck  — D~3Y>ty(y2ck)\ 
is  augmented  but  never  above  the  value  of  the  error  |yicfc  ~  D~3Y,ty(yick)\  which  is 
evident  in  Figure  4.35.  The  optimal  valve  commands  U\  and  u2  are  plotted  in  Fig- 
ure 4.36  and  Figure  4.37  for  A  =  2  and  A  =  3,  respectively.  For  A  =  2,  the  signal  U\ 
is  constrained  throughout  most  of  the  simulation  while  for  A  =  3  the  signal  W2,  which 
influences  y\  in  three  steps,  is  optimally  used  to  reduce  the  error  \y\ck  —  D~3Y,ty(yick)\- 


CHAPTER  5 
CONCLUSION 


5.1     Summary 

This  dissertation  presents  new  results  for  the  problem  of  disturbance  atten- 
uation for  nonlinear  control  systems.  Of  primary  significance  is  the  derivation  of  a 
performance  bound  that  provides  a  measure  of  the  ability  to  match  a  set  by  subset 
of  the  image  of  a  system.  Our  bound  arose  from  the  desire  of  providing  an  estimate 
of  the  minimal  effect  of  the  disturbance  signal  d  on  the  output  signal  z  of  Figure  1.1. 

It  turns  out  that  the  calculation  of  the  performance  bound  is  rather  labori- 
ous because  its  definition  (Definition  3.4.1)  uses  ImS.  In  practice,  we  will  find  an 
estimate  of  the  performance  bound  that  uses  an  approximate  right  inverse  system 
in  its  calculation  instead  of  ImE.  The  estimate  of  the  performance  bound  involves 
an  optimization  problem  which  relates  to  finding  a  global  maximum  of  a  non-convex 
function.  Nonlinear  programming  is  used  to  calculate  the  estimate  of  the  performance 
bound. 

Three  systems  of  practical  origin  were  selected  for  the  analysis  of  the  per- 
formance bound.  In  the  PM  stepper  motor  example,  a  closed  form  estimate  of  the 
performance  bound  was  found  which  was  verified  by  simulation.  In  the  aerodynamic 
model  example,  a  nonlinear  programming  algorithm  was  determined  to  compute  the 
estimate  of  the  performance  bound.  The  multivariable  process  model  example  pro- 
vided us  insight  into  how  different  design  integers  can  influence  the  estimate  of  the 
performance  bound. 
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5.2     Future  Directions 

The  results  of  this  dissertation  serve  to  enhance  many  of  the  tools  used  to 
obtain  the  optimal  controller  for  nonlinear  control  systems.  In  the  process  of  com- 
pleting one  particular  objective,  other  promising  areas  frequently  appear  as  directions 
for  additional  research.  One  of  the  more  intriguing  of  such  areas  emerging  from  this 
work  is  the  derivation  of  a  bound  that  handles  constant  bias  disturbances. 

Referring  to  Figure  1.1,  and  using  the  parameterization  (1.3)  with  the  assump- 
tion that  the  system  E  is  stable  (i.e.,  P  —  E  and  Q  =  /),  we  get 

Y,c{v,d)  =  d+Y.<l>{v,d).  (5.1) 

The  stable  and  causal  system  <J>  is  selected  such  that 

(j>(v,d)  =  ${Zv-d)  (5.2) 

where  $  is  an  approximate  right  inverse  system.  All  that  is  needed  to  define  a 
performance  bound  is  an  optimization  set.  In  this  case,  the  optimization  set  is  the 
image  of  [Eu  —  d]  where  the  domain  of  [Et>  —  d]  is  a  bounded  set.  The  performance 
bound  is  given  by 

-P(E)  =  gt  (im  [Ev  -  d] ,  £rAImE)  .  (5.3) 

The  bound  of  (5.3)  is  a  theoretical  estimate  of  the  minimal  effect  of  the  disturbance 
d  on  the  output  z. 

Let's  use  A  =  0  initially  in  the  bound  optimization  problem.  Note  that  there 
is  another  A  associated  with  the  approximate  inverse  optimization  problem.  Both  are 
equal  for  most  applications.  In  this  application  the  design  integer  A  will  be  associated 
with  an  approximate  inverse  and  A  >  0.  Now  we  need  to  modify  our  definition  of  the 
norm  pt  for  the  bound.  For  a  sequence  u  6  S{Wn)  and  for  a  j  >  0,  the  norm  /t»ej  is 
defined 

peJ(u)  =%up  (1  +  eP  |ut|,     6>0.  (5.4) 

i>j 
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The  class  of  disturbances  is  a  constant  bias;  i.e.,  the  magnitude  of  the  disturbance  is 
unknown,  but  it  is  constant  (i.e.,  a  step  disturbance).  Let  us  define  our  performance 
bound  for  A  =  0,  denoted  by  7£(E),  such  that 

7£(E)  =  inf  sup  pt<\  (d  +  E<f)(v,d)  —  Ei>)  (5.5) 

where  £t>  is  our  desired  response.  The  supremum  is  over  the  input  signal  space  and 
the  infimum  is  over  the  stable  and  causal  systems  4>.  The  stable  and  causal  system 
4>  can  be  chosen  as 

<t>(v,d)  =  $(D-xEv-d)  (5.6) 

where  the  design  integer  A  >  0  of  an  approximate  right  inverse  system  ^  is  the  least 
latency  of  the  system  E. 

It  will  be  shown  that  the  bound  7£(E)  is  well  suited  for  a  constant  bias  dis- 
turbance d.  If  the  system  ty  were  a  right  inverse  of  the  system  E,  then  E^  =  D  I. 
Let  us  substitute  (5.6)  into  (5.5),  we  get 

K(E)    =  sup/>e,A  (d  +  Etf  (D~xEv  -d)  -  Ei>) 

=  sup /?e,A  (d  +  DXI  (D~xZv  -  d)  -  £i>) 

=  sup  /3e,A  (d  +  (Ev  -  Dxd)  -  Ev) 

=  sup  pet\  (d  —  D  d) 

11(E)    =  0 

where  the  supremum  is  taken  over  the  input  signal  space. 

The  bound  72.(E)  can  be  calculated  using  techniques  of  Chapter  4.  The  op- 
timization set  is  the  image  of  [D~xIlv  —  d]  where  the  domain  of  [D~xHv  —  d\  is  a 
bounded  set.  For  the  PM  stepper  motor  example,  the  estimate  of  the  performance 
bound  'P(S)  of  (4.26)  turns  out  to  be  the  performance  bound  7£(£)  because  the 
performance  bound  for  this  problem  is  independent  of  the  optimization  set.  Thus 

K{Y)  =  ^y^-^t2  (5-7) 
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for  the  PM  stepper  motor  example. 

Throughout  this  dissertation,  no  new  results  were  found  for  the  design  of  the 
equivalent  controller  C  of  Figure  1.1  because  the  results  from  [22]  were  sufficient.  In 
this  emerging  area  of  of  deriving  a  performance  bound  7£(E)  to  handle  a  constant 
bias,  attention  should  be  given  to  the  controllers  that  implement  the  new  control  laws. 
There  will  be  other  performance  requirements  in  addition  to  the  bound  requirement, 
and  so  the  performance  index  will  be  augmented  to  include  additional  variables. 
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